{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import pymc4 as pm\n",
    "import arviz as az\n",
    "\n",
    "import tensorflow_probability as tfp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<pymc4.random_variables.continuous.HalfNormal at 0x1c3d18dc88>,\n",
       " <pymc4.random_variables.continuous.Normal at 0x1c3d18d940>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@pm.model(auto_name=True)\n",
    "def t_test():\n",
    "    sigma = pm.HalfNormal(1)\n",
    "    mu = pm.Normal(0, sigma)\n",
    "\n",
    "model = t_test.configure()\n",
    "\n",
    "model._forward_context.vars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: id=109, shape=(), dtype=float32, numpy=-3.5687468>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model._forward_context.vars[0]._distribution.sample()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'sigma': <tf.Tensor: id=136, shape=(), dtype=float32, numpy=1.1835041>,\n",
       " 'mu': <tf.Tensor: id=158, shape=(), dtype=float32, numpy=-1.2197603>}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.forward_sample()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tensorflow.python.eager.def_function.Function at 0x1c3d18d588>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.function(model.make_log_prob_function())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: id=290, shape=(), dtype=float32, numpy=-4.5149107>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func = model.make_log_prob_function()\n",
    "\n",
    "sd = -tf.ones((1,))\n",
    "mu = tf.ones((1,))\n",
    "\n",
    "func(sd, mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample(model, nsteps=200, burnin=100, leapfrog_steps=10, compile=True):\n",
    "    # Since HMC operates over unconstrained space, we need to transform the\n",
    "    # samples so they live in real-space.\n",
    "    random_variables = model._forward_context.vars\n",
    "    unconstraining_bijectors = []\n",
    "    inits = []\n",
    "    for rv in random_variables:\n",
    "        inits.append(tf.ones(rv.sample().shape))\n",
    "        unconstraining_bijectors.append(tfp.bijectors.Identity())\n",
    "\n",
    "    unnormalized_posterior_log_prob = model.make_log_prob_function()\n",
    "    \n",
    "    if compile:\n",
    "        # compile logp to speed things up\n",
    "        unnormalized_posterior_log_prob = tf.function(\n",
    "            unnormalized_posterior_log_prob)\n",
    "\n",
    "    # Initialize the step_size. (It will be automatically adapted.)\n",
    "#     step_size = tf.get_variable(\n",
    "#         name='step_size',\n",
    "#         initializer=tf.constant(0.5, dtype=tf.float32),\n",
    "#         trainable=False,\n",
    "#         use_resource=True\n",
    "#     )\n",
    "\n",
    "    sample_chain = tfp.mcmc.sample_chain\n",
    "    \n",
    "    # Defining the HMC\n",
    "    hmc = tfp.mcmc.TransformedTransitionKernel(\n",
    "        inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(\n",
    "            target_log_prob_fn=unnormalized_posterior_log_prob,\n",
    "            num_leapfrog_steps=leapfrog_steps,\n",
    "            step_size=.5,\n",
    "            step_size_update_fn=tfp.mcmc.make_simple_step_size_update_policy(100),\n",
    "            state_gradients_are_stopped=False),\n",
    "        bijector=unconstraining_bijectors)\n",
    "\n",
    "    # Sampling from the chain.\n",
    "    posterior_samples_tensor, kernel_results = sample_chain(\n",
    "        num_results=nsteps,\n",
    "        num_burnin_steps=burnin,\n",
    "        current_state=inits,\n",
    "        kernel=hmc)\n",
    "\n",
    "    # Initialize any created variables.\n",
    "    init_g = tf.global_variables_initializer()\n",
    "    init_l = tf.local_variables_initializer()\n",
    "    \n",
    "    evaluate(init_g)\n",
    "    evaluate(init_l)\n",
    "    \n",
    "    *posterior_samples, kernel_results_ = evaluate([\n",
    "        *posterior_samples_tensor,\n",
    "        kernel_results,\n",
    "    ])\n",
    "    \n",
    "    trace = {rv.name: arr for rv, arr in zip(random_variables, posterior_samples)}\n",
    "    \n",
    "    return az.dict_to_dataset(trace), kernel_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/sample.py:335: UserWarning: Tracing all kernel results by default is deprecated. Set the `trace_fn` argument to None (the future default value) or an explicit callback that traces the values you are interested in.\n",
      "  warnings.warn(\"Tracing all kernel results by default is deprecated. Set \"\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'float' object has no attribute 'assign_add'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m<ipython-input-13-bb5140909246>\u001b[0m in \u001b[0;36msample\u001b[0;34m(model, nsteps, burnin, leapfrog_steps, compile)\u001b[0m\n\u001b[1;32m     41\u001b[0m         \u001b[0mnum_burnin_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mburnin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     42\u001b[0m         \u001b[0mcurrent_state\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minits\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 43\u001b[0;31m         kernel=hmc)\n\u001b[0m\u001b[1;32m     44\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     45\u001b[0m     \u001b[0;31m# Initialize any created variables.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/sample.py\u001b[0m in \u001b[0;36msample_chain\u001b[0;34m(num_results, current_state, previous_kernel_results, kernel, num_burnin_steps, num_steps_between_results, trace_fn, return_final_kernel_results, parallel_iterations, name)\u001b[0m\n\u001b[1;32m    359\u001b[0m                                             trace_fn(*state_and_results)),\n\u001b[1;32m    360\u001b[0m         \u001b[0;31m# pylint: enable=g-long-lambda\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 361\u001b[0;31m         parallel_iterations=parallel_iterations)\n\u001b[0m\u001b[1;32m    362\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    363\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mreturn_final_kernel_results\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36mtrace_scan\u001b[0;34m(loop_fn, initial_state, elems, trace_fn, parallel_iterations, name)\u001b[0m\n\u001b[1;32m    367\u001b[0m         \u001b[0mbody\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m_body\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    368\u001b[0m         \u001b[0mloop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_state\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrace_arrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 369\u001b[0;31m         parallel_iterations=parallel_iterations)\n\u001b[0m\u001b[1;32m    370\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    371\u001b[0m     \u001b[0mstacked_trace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmap_structure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrace_arrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mwhile_loop_v2\u001b[0;34m(cond, body, loop_vars, shape_invariants, parallel_iterations, back_prop, swap_memory, maximum_iterations, name)\u001b[0m\n\u001b[1;32m   3219\u001b[0m       \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3220\u001b[0m       \u001b[0mmaximum_iterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaximum_iterations\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3221\u001b[0;31m       return_same_structure=True)\n\u001b[0m\u001b[1;32m   3222\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3223\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mwhile_loop\u001b[0;34m(cond, body, loop_vars, shape_invariants, parallel_iterations, back_prop, swap_memory, name, maximum_iterations, return_same_structure)\u001b[0m\n\u001b[1;32m   3445\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3446\u001b[0m       \u001b[0;32mwhile\u001b[0m \u001b[0mcond\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mloop_vars\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3447\u001b[0;31m         \u001b[0mloop_vars\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbody\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mloop_vars\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3448\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mtry_to_pack\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mloop_vars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_basetuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3449\u001b[0m           \u001b[0mpacked\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36m_body\u001b[0;34m(i, state, trace_arrays)\u001b[0m\n\u001b[1;32m    356\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    357\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_body\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrace_arrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 358\u001b[0;31m       \u001b[0mstate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mloop_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0melems_array\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    359\u001b[0m       trace_arrays = tf.nest.pack_sequence_as(trace_arrays, [\n\u001b[1;32m    360\u001b[0m           a.write(i, v) for a, v in zip(\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/sample.py\u001b[0m in \u001b[0;36m_trace_scan_fn\u001b[0;34m(state_and_results, num_steps)\u001b[0m\n\u001b[1;32m    343\u001b[0m           \u001b[0mbody_fn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mone_step\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m           \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstate_and_results\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m           parallel_iterations=parallel_iterations)\n\u001b[0m\u001b[1;32m    346\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0mnext_state\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcurrent_kernel_results\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36msmart_for_loop\u001b[0;34m(loop_num_iter, body_fn, initial_loop_vars, parallel_iterations, name)\u001b[0m\n\u001b[1;32m    283\u001b[0m           \u001b[0mbody\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbody_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    284\u001b[0m           \u001b[0mloop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 285\u001b[0;31m           \u001b[0mparallel_iterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mparallel_iterations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    286\u001b[0m       )[1:]\n\u001b[1;32m    287\u001b[0m     \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mwhile_loop_v2\u001b[0;34m(cond, body, loop_vars, shape_invariants, parallel_iterations, back_prop, swap_memory, maximum_iterations, name)\u001b[0m\n\u001b[1;32m   3219\u001b[0m       \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3220\u001b[0m       \u001b[0mmaximum_iterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaximum_iterations\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3221\u001b[0;31m       return_same_structure=True)\n\u001b[0m\u001b[1;32m   3222\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3223\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mwhile_loop\u001b[0;34m(cond, body, loop_vars, shape_invariants, parallel_iterations, back_prop, swap_memory, name, maximum_iterations, return_same_structure)\u001b[0m\n\u001b[1;32m   3445\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3446\u001b[0m       \u001b[0;32mwhile\u001b[0m \u001b[0mcond\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mloop_vars\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3447\u001b[0;31m         \u001b[0mloop_vars\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbody\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mloop_vars\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3448\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mtry_to_pack\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mloop_vars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_basetuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3449\u001b[0m           \u001b[0mpacked\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(i, *args)\u001b[0m\n\u001b[1;32m    281\u001b[0m       return tf.while_loop(\n\u001b[1;32m    282\u001b[0m           \u001b[0mcond\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mloop_num_iter\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 283\u001b[0;31m           \u001b[0mbody\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbody_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    284\u001b[0m           \u001b[0mloop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    285\u001b[0m           \u001b[0mparallel_iterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mparallel_iterations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/transformed_kernel.py\u001b[0m in \u001b[0;36mone_step\u001b[0;34m(self, current_state, previous_kernel_results)\u001b[0m\n\u001b[1;32m    256\u001b[0m       transformed_next_state, kernel_results = self._inner_kernel.one_step(\n\u001b[1;32m    257\u001b[0m           \u001b[0mprevious_kernel_results\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransformed_state\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 258\u001b[0;31m           previous_kernel_results.inner_results)\n\u001b[0m\u001b[1;32m    259\u001b[0m       transformed_next_state_parts = (\n\u001b[1;32m    260\u001b[0m           \u001b[0mtransformed_next_state\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/hmc.py\u001b[0m in \u001b[0;36mone_step\u001b[0;34m(self, current_state, previous_kernel_results)\u001b[0m\n\u001b[1;32m    535\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep_size_update_fn\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    536\u001b[0m         step_size_assign = self.step_size_update_fn(  # pylint: disable=not-callable\n\u001b[0;32m--> 537\u001b[0;31m             self.step_size, kernel_results)\n\u001b[0m\u001b[1;32m    538\u001b[0m         kernel_results = kernel_results._replace(\n\u001b[1;32m    539\u001b[0m             extra=HamiltonianMonteCarloExtraKernelResults(\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/hmc.py\u001b[0m in \u001b[0;36mstep_size_simple_update_fn\u001b[0;34m(step_size_var, kernel_results)\u001b[0m\n\u001b[1;32m    158\u001b[0m             \u001b[0mpred\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstep_counter\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mnum_adaptation_steps\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    159\u001b[0m             \u001b[0mtrue_fn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mbuild_assign_op\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 160\u001b[0;31m             false_fn=lambda: step_size_var)\n\u001b[0m\u001b[1;32m    161\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    162\u001b[0m   \u001b[0;32mreturn\u001b[0m \u001b[0mstep_size_simple_update_fn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mcond_for_tf_v2\u001b[0;34m(pred, true_fn, false_fn, name)\u001b[0m\n\u001b[1;32m   2145\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2146\u001b[0m   \"\"\"\n\u001b[0;32m-> 2147\u001b[0;31m   \u001b[0;32mreturn\u001b[0m \u001b[0mcond\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpred\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrue_fn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtrue_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfalse_fn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfalse_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstrict\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2148\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2149\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/util/deprecation.py\u001b[0m in \u001b[0;36mnew_func\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    505\u001b[0m                 \u001b[0;34m'in a future version'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'after %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mdate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    506\u001b[0m                 instructions)\n\u001b[0;32m--> 507\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    508\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    509\u001b[0m     doc = _add_deprecated_arg_notice_to_docstring(\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mcond\u001b[0;34m(pred, true_fn, false_fn, strict, name, fn1, fn2)\u001b[0m\n\u001b[1;32m   1944\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mcontext\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexecuting_eagerly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1945\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0mpred\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1946\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_UnpackIfSingleton\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrue_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1947\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0m_UnpackIfSingleton\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfalse_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1948\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/tensorflow_probability/python/mcmc/hmc.py\u001b[0m in \u001b[0;36mbuild_assign_op\u001b[0;34m()\u001b[0m\n\u001b[1;32m    148\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0mss\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mstep_size_var\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    149\u001b[0m         ]\n\u001b[0;32m--> 150\u001b[0;31m       return step_size_var.assign_add(\n\u001b[0m\u001b[1;32m    151\u001b[0m           step_size_var * tf.cast(adjustment, step_size_var.dtype))\n\u001b[1;32m    152\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'float' object has no attribute 'assign_add'"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "trace, sampler_stats = sample(model, nsteps=1000, compile=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (chain: 1, draw: 1000)\n",
       "Coordinates:\n",
       "  * chain    (chain) int64 0\n",
       "  * draw     (draw) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...\n",
       "Data variables:\n",
       "    sigma    (chain, draw) float32 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ...\n",
       "    mu       (chain, draw) float32 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ...\n",
       "Attributes:\n",
       "    created_at:  2019-03-11T10:10:16.355586"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:2507: RuntimeWarning: invalid value encountered in true_divide\n",
      "  pk = 1.0*pk / np.sum(pk, axis=0)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "cannot convert float NaN to integer",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-7c3f6fe9d699>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0maz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot_trace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/arviz/plots/traceplot.py\u001b[0m in \u001b[0;36mplot_trace\u001b[0;34m(data, var_names, coords, figsize, textsize, lines, combined, kde_kwargs, hist_kwargs, trace_kwargs)\u001b[0m\n\u001b[1;32m    117\u001b[0m                 \u001b[0m_histplot_op\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhist_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    118\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 119\u001b[0;31m                 \u001b[0mplot_kde\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtextsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtextsize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkde_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    120\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    121\u001b[0m         \u001b[0maxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_yticks\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/arviz/plots/kdeplot.py\u001b[0m in \u001b[0;36mplot_kde\u001b[0;34m(values, values2, cumulative, rug, label, bw, rotated, contour, fill_last, figsize, textsize, plot_kwargs, fill_kwargs, rug_kwargs, contour_kwargs, ax)\u001b[0m\n\u001b[1;32m    136\u001b[0m         \u001b[0mrug_kwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msetdefault\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'markersize'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mmarkersize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    137\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 138\u001b[0;31m         \u001b[0mdensity\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlower\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupper\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_fast_kde\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcumulative\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    139\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlower\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupper\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdensity\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    140\u001b[0m         \u001b[0mfill_func\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfill_between\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/arviz/plots/kdeplot.py\u001b[0m in \u001b[0;36m_fast_kde\u001b[0;34m(x, cumulative, bw)\u001b[0m\n\u001b[1;32m    210\u001b[0m     \u001b[0mstd_x\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mentropy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mxmin\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mbw\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    211\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 212\u001b[0;31m     \u001b[0mn_bins\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen_x\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mstd_x\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m200\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    213\u001b[0m     \u001b[0mgrid\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhistogram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbins\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mn_bins\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    214\u001b[0m     \u001b[0md_x\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mxmax\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mxmin\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mn_bins\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: cannot convert float NaN to integer"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAD8CAYAAACxZPjGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH4dJREFUeJzt3X+QXXd53/H3JxK2g0OwQEqHWjayiQA70LFhx5DSCSRgI9yORRsmkSnFUCcaCCZtaKdjhpmQmj9CSAoMUwesBJUf09iAmyEbxozHAXvoAAKtimNsEWMhCN7IjUVs3E5NbGw//eMctderXe05u/eH7u77NXNn7/me77n7fHXuPvfRud9zTqoKSZIkSd39xKQDkCRJkqaNRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1NOyRXSSvUnuT3LnEuuT5ENJDiW5I8mLBtZdkeSe9nHFMAOXJEmSJqXLkeiPATtOsP41wPb2sRv4MECSZwDvBl4CXAS8O8mm1QQrSZIknQyWLaKr6kvAAyfoshP4RDX2AWckeRbwauCWqnqgqh4EbuHExbgkSZI0FTYO4TXOBO4dWJ5v25ZqP06S3TRHsTn99NNf/PznP38IYUnS+B04cOAHVbVl0nGM0+bNm2vbtm2TDkOSeltNzh5GEZ1F2uoE7cc3Vu0B9gDMzMzU3NzcEMKSpPFL8teTjmHctm3bhnlb0jRaTc4extU55oGzBpa3AkdO0C5JGqFVnhD+eJLb28fs+KKWpOkyjCJ6Fnhjm5RfCjxUVfcBNwOXJNnUnlB4SdsmSRqtj7GCE8JbP6qqC9rHZaMLUZKm27LTOZJcD7wC2JxknuaKG08BqKqPADcBlwKHgIeBN7frHkjyHmB/+1LXVNWJTlCUJA1BVX0pybYTdPl/J4QD+5KckeRZ7QEQSVIHyxbRVXX5MusLeNsS6/YCe1cWmiRpRJY68fs+4LQkc8BjwHur6rOLvcDgCeFnn332aKOVpJOQdyyUpPXnRCd+n11VM8DrgQ8mec5iL1BVe6pqpqpmtmxZVxcjkSTAIlqS1qMlT/yuqmM/DwO3AReOOzhJmgYW0ZK0/ix6Qnh7IvipAEk2Ay8DDk4yUEk6WQ3jOtGSpJPISk8IB84DrkvyBM1BlvdWlUW0JC3CIlqS1piVnhBeVV8BXjiquCRpLXE6hyRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUU6ciOsmOJHcnOZTk6kXWfyDJ7e3j20l+OLDu8YF1s8MMXpIkSZqEjct1SLIBuBa4GJgH9ieZraqDx/pU1W8N9H87cOHAS/yoqi4YXsiSJEnSZHU5En0RcKiqDlfVo8ANwM4T9L8cuH4YwUmSJEknoy5F9JnAvQPL823bcZI8GzgH+OJA82lJ5pLsS/LaJbbb3faZO3r0aMfQJUmSpMnoUkRnkbZaou8u4Maqenyg7eyqmgFeD3wwyXOOe7GqPVU1U1UzW7Zs6RCSJEmSNDldiuh54KyB5a3AkSX67mLBVI6qOtL+PAzcxpPnS0uSJElTp0sRvR/YnuScJKfQFMrHXWUjyfOATcBXB9o2JTm1fb4ZeBlwcOG2kiRJ0jRZ9uocVfVYkquAm4ENwN6quivJNcBcVR0rqC8Hbqiqwake5wHXJXmCpmB/7+BVPSRJkqRptGwRDVBVNwE3LWj77QXLv7PIdl8BXriK+CRJkqSTjncslCRJknqyiJYkSZJ6soiWJEmSerKIliRJknqyiJYkSZJ6soiWJEmSerKIliRJknqyiJYkSZJ6soiWJEmSerKIliRJknqyiJYkSZJ6soiWJEmSerKIliRJknqyiJYkSZJ6soiWJEmSerKIliRJknrqVEQn2ZHk7iSHkly9yPo3JTma5Pb28WsD665Ick/7uGKYwUuSjpdkb5L7k9y5xPok+VCb0+9I8qKBdeZsSepg43IdkmwArgUuBuaB/Ulmq+rggq6fqqqrFmz7DODdwAxQwIF22weHEr0kaTEfA/4z8Ikl1r8G2N4+XgJ8GHiJOVuSulu2iAYuAg5V1WGAJDcAO4GFRfRiXg3cUlUPtNveAuwArl9ZuJKk5VTVl5JsO0GXncAnqqqAfUnOSPIs4BWMKWff+TcP8dCPfjzsl5W0Tj39J5/CC858+lh/Z5fpHGcC9w4sz7dtC/1y+7XgjUnO6rNtkt1J5pLMHT16tGPokqQVWio3d8335m1J616XI9FZpK0WLP85cH1VPZLkLcDHgV/quC1VtQfYAzAzM3PceknSUC2VmzvlbFh93h73ESNJGrYuR6LngbMGlrcCRwY7VNXfVdUj7eIfAS/uuq0kaeyWys3mbEnqqEsRvR/YnuScJKcAu4DZwQ7tXLpjLgO+1T6/GbgkyaYkm4BL2jZJ0uTMAm9sr9LxUuChqroPc7YkdbbsdI6qeizJVTSJdAOwt6ruSnINMFdVs8BvJrkMeAx4AHhTu+0DSd5DU4gDXHPshBVJ0mgkuZ7mJMHNSeZprrjxFICq+ghwE3ApcAh4GHhzu86cLUkdpTk5++QxMzNTc3Nzkw5DklYkyYGqmpl0HONk3pY0rVaTs71joSRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1FOnIjrJjiR3JzmU5OpF1r8jycEkdyT5QpJnD6x7PMnt7WN2mMFLkiRJk7BxuQ5JNgDXAhcD88D+JLNVdXCg2zeAmap6OMlbgfcBv9qu+1FVXTDkuCVJkqSJ6XIk+iLgUFUdrqpHgRuAnYMdqurWqnq4XdwHbB1umJIkSdLJo0sRfSZw78DyfNu2lCuBzw8sn5ZkLsm+JK9dbIMku9s+c0ePHu0QkiRJkjQ5y07nALJIWy3aMXkDMAO8fKD57Ko6kuRc4ItJvllV33nSi1XtAfYAzMzMLPrakiRJ0smiy5HoeeCsgeWtwJGFnZK8CngXcFlVPXKsvaqOtD8PA7cBF64iXkmSJGniuhTR+4HtSc5JcgqwC3jSVTaSXAhcR1NA3z/QvinJqe3zzcDLgMETEiVJkqSps+x0jqp6LMlVwM3ABmBvVd2V5Bpgrqpmgd8Hfgr4TBKA71fVZcB5wHVJnqAp2N+74KoekiRJ0tTpMieaqroJuGlB228PPH/VEtt9BXjhagKUJEmSTjbesVCSJEnqySJakiRJ6skiWpIkSerJIlqSJEnqySJakiRJ6skiWpIkSerJIlqSJEnqySJakiRJ6skiWpIkSerJIlqSJEnqySJakiRJ6skiWpIkSerJIlqS1pgkO5LcneRQkqsXWf/sJF9IckeS25JsHVj3eJLb28fseCOXpOmxcdIBSJKGJ8kG4FrgYmAe2J9ktqoODnT7A+ATVfXxJL8E/C7wr9p1P6qqC8YatCRNIY9ES9LachFwqKoOV9WjwA3AzgV9zge+0D6/dZH1kqRlWERL0tpyJnDvwPJ82zboL4Ffbp//c+BpSZ7ZLp+WZC7JviSvXeqXJNnd9ps7evTosGKXpKnRqYjuML/u1CSfatd/Lcm2gXXvbNvvTvLq4YUuSVpEFmmrBcv/Hnh5km8ALwf+BnisXXd2Vc0Arwc+mOQ5i/2SqtpTVTNVNbNly5YhhS5J02PZInpgft1raL4CvDzJ+Qu6XQk8WFU/C3wA+L122/OBXcDPATuAP2xfT5I0GvPAWQPLW4Ejgx2q6khV/YuquhB4V9v20LF17c/DwG3AhWOIWZKmTpcj0V3m1+0EPt4+vxF4ZZK07TdU1SNV9V3gUPt6kqTR2A9sT3JOklNoDmQ86SobSTYnOZb/3wnsbds3JTn1WB/gZcDgCYmSpFaXq3MsNr/uJUv1qarHkjwEPLNt37dg24Vz80iyG9jdLj6S5M5O0a8dm4EfTDqIMXPM68N6HPPzJvnL2xx8FXAzsAHYW1V3JbkGmKuqWeAVwO8mKeBLwNvazc8DrkvyBM1BlvcuuKrHog4cOPCDJH+9gnDX8vvDsU2vtTw+x3a8Z6/0F3YporvMr1uqT5dtqao9wB6AJHPtfLx1wzGvD455fUgyN+kYquom4KYFbb898PxGmm8NF273FeCFK/h9K5oUvZbfH45teq3l8Tm24eoynWPZ+XWDfZJsBJ4OPNBxW0mSJGmqdCmil51f1y5f0T5/HfDFqqq2fVd79Y5zgO3A14cTuiRJkjQZy07n6Di/7qPAJ5McojkCvavd9q4kn6Y5MeUx4G1V9fgyv3LPyocztRzz+uCY14f1OOaVWsv/Vo5teq3l8Tm2IUpzwFiSJElSV96xUJIkSerJIlqSJEnqaWJF9GpuJT6tOoz5HUkOJrkjyReSrPjahSeL5cY80O91SSrJ1F96p8uYk/xKu6/vSvIn445x2Dq8t89OcmuSb7Tv70snEeewJNmb5P6lrmmfxofaf487krxo3DGezLrmhZNVkrPa9/O32r/hf9O2PyPJLUnuaX9uatun7v2QZEP79/q5dvmc9rP4nvaz+ZS2feo+q5OckeTGJH/V7sOfXyv7Lslvte/JO5Ncn+S0ad53i+XaleyrJFe0/e9JcsViv2tFqmrsD5oTFL8DnAucAvwlcP6CPr8BfKR9vgv41CRiHfOYfxF4avv8rethzG2/p9Hc8GEfMDPpuMewn7cD3wA2tcs/M+m4xzDmPcBb2+fnA9+bdNyrHPMvAC8C7lxi/aXA52mulf9S4GuTjvlkeXTNCyfzA3gW8KL2+dOAb7fv6/cBV7ftVwO/N63vB+AdwJ8An2uXPw3sap9/ZODveeo+q2nusPxr7fNTgDPWwr6juZndd4GfHNhnb5rmfbdYru27r4BnAIfbn5va55uGEd+kjkSv5lbi02rZMVfVrVX1cLu4j+a62tOsy34GeA/NH8XfjzO4Eeky5l8Hrq2qBwGq6v4xxzhsXcZcwE+3z5/OlF8vvqq+RHMloqXsBD5RjX3AGUmeNZ7oTnpd88JJq6ruq6r/0T7/38C3aAqYwc+tjwOvbZ9P1fshyVbgnwJ/3C4H+CX+/w16Fo5taj6rk/w0TWH2UYCqerSqfsga2Xc0V137yTT37HgqcB9TvO+WyLV999WrgVuq6oH2c/cWYMcw4ptUEb3YrcQX3g78SbcSB47dSnxadRnzoCtp/kc1zZYdc5ILgbOq6nPjDGyEuuzn5wLPTfLlJPuSDOWPeYK6jPl3gDckmae5k97bxxPaxPT9e19P1tS/TfsV+IXA14B/UFX3QVNoAz/Tdpu2MX8Q+A/AE+3yM4Eftp/F8OT4p+2z+lzgKPBf2ukqf5zkdNbAvquqvwH+APg+TfH8EHCAtbPvjum7r0a2DydVRK/mVuLTqvN4krwBmAF+f6QRjd4Jx5zkJ4APAP9ubBGNXpf9vJFmSscrgMuBP05yxojjGqUuY74c+FhVbaX5yu2T7f5fq9Za/hqmNfNvk+SngP8G/Nuq+l8n6rpI20k55iT/DLi/qg4MNi/StTqsOxltpJke8OGquhD4PzRTApYyNeNr5wbvBM4B/iFwOvCaRbpO675bzlLjGdk4J/UhtppbiU+rTrdAT/Iq4F3AZVX1yJhiG5Xlxvw04AXAbUm+RzOHaTbTfXJh1/f2n1XVj6vqu8DdNEX1tOoy5itp5uVRVV8FTgM2jyW6yej0975OrYl/myRPoSmg/2tV/Wnb/LfHvupvfx6bqjVNY34ZcFmbk2+gmQrwQZqvxo/doG0w/mn7rJ4H5qvqa+3yjTRF9VrYd68CvltVR6vqx8CfAv+YtbPvjum7r0a2DydVRK/mVuLTatkxt1MbrqMpoKd9niwsM+aqeqiqNlfVtqraRjMP/LKqmptMuEPR5b39WZqTSEmymWZ6x+GxRjlcXcb8feCVAEnOoymij441yvGaBd7Yni3+UuChY18/qtP75aTWzhv9KPCtqnr/wKrBz60rgD8baJ+K90NVvbOqtrY5eRfNZ++/BG6l+SyG48c2NZ/VVfU/gXuTPK9teiXNXZWnft/R5NmXJnlq+x49NrY1se8G9N1XNwOXJNnUHq2/pG1bvWGcnbiSB81Xut+mOUv7XW3bNTRFFDQfsp8BDgFfB86dVKxjHPNfAH8L3N4+Zicd86jHvKDvbUz51Tk67ucA76dJbt+kPWt6mh8dxnw+8GWaKzHcDlwy6ZhXOd7raeYc/pjmKMeVwFuAtwzs42vbf49vroX39ajfL9P0AP4JzdfBdwzk60tp5pN+Abin/fmMaX4/0Ew5O3Z1jnPbz+JD7WfzqW371H1WAxcAc+3++yzNFRvWxL4D/iPwV8CdwCeBU6d53y2Ra3vvK+Bft+M8BLx5WPF5229JkiSpp2Wncyx2oesF68d/cWtJ0pLM25I0el3mRH+ME19P7zU0J0VtB3YDH4bmjjLAu4GX0FwX9N3tXBRJ0mh9DPO2JI3UskV0rfymAiO7uLUkaWnmbUkavY3Ld1nWqi9unWQ3zdEQTj/99Bc///nPH0JYkjR+Bw4c+EFVbZl0HMswb0sSq8vZwyiiV31x66raA+wBmJmZqbm5ab7CmaT1LMlfTzqGDszbksTqcvYwrhM99otbS5JWxbwtSas0jCJ6/Be3liSthnlbklZp2ekcSa6nueD65iTzNGduPwWgqj4C3ERzkflDwMPAm9t1DyR5D83dqQCuqappuJ2kJE0187Ykjd6yRXRVXb7M+gLetsS6vcDelYUmSVoJ87Ykjd4wpnNIkiRJ64pFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1FOnIjrJjiR3JzmU5OpF1n8gye3t49tJfjiw7vGBdbPDDF6SdDxztiSN3sblOiTZAFwLXAzMA/uTzFbVwWN9quq3Bvq/Hbhw4CV+VFUXDC9kSdJSzNmSNB5djkRfBByqqsNV9ShwA7DzBP0vB64fRnCSpN7M2ZI0Bl2K6DOBeweW59u24yR5NnAO8MWB5tOSzCXZl+S1S2y3u+0zd/To0Y6hS5IWMfKc3W5r3pa0rnUporNIWy3RdxdwY1U9PtB2dlXNAK8HPpjkOce9WNWeqpqpqpktW7Z0CEmStISR52wwb0tSlyJ6HjhrYHkrcGSJvrtY8LVgVR1pfx4GbuPJc+8kScNlzpakMehSRO8Htic5J8kpNEn3uDO2kzwP2AR8daBtU5JT2+ebgZcBBxduK0kaGnO2JI3BslfnqKrHklwF3AxsAPZW1V1JrgHmqupYcr4cuKGqBr82PA+4LskTNAX7ewfPEJckDZc5W5LGI0/On5M3MzNTc3Nzkw5DklYkyYF2TvG6Yd6WNK1Wk7O9Y6EkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktSTRbQkSZLUk0W0JEmS1JNFtCRJktRTpyI6yY4kdyc5lOTqRda/KcnRJLe3j18bWHdFknvaxxXDDF6SdDxztiSN3sblOiTZAFwLXAzMA/uTzFbVwQVdP1VVVy3Y9hnAu4EZoIAD7bYPDiV6SdKTmLMlaTy6HIm+CDhUVYer6lHgBmBnx9d/NXBLVT3QJuFbgB0rC1WS1IE5W5LGoEsRfSZw78DyfNu20C8nuSPJjUnO6rNtkt1J5pLMHT16tGPokqRFjDxng3lbkroU0VmkrRYs/zmwrar+EfAXwMd7bEtV7amqmaqa2bJlS4eQJElLGHnOBvO2JHUpoueBswaWtwJHBjtU1d9V1SPt4h8BL+66rSRpqMzZkjQGXYro/cD2JOckOQXYBcwOdkjyrIHFy4Bvtc9vBi5JsinJJuCStk2SNBrmbEkag2WvzlFVjyW5iiaRbgD2VtVdSa4B5qpqFvjNJJcBjwEPAG9qt30gyXtokjrANVX1wAjGIUnCnC1J45KqRae7TczMzEzNzc1NOgxJWpEkB6pqZtJxjJN5W9K0Wk3O9o6FkiRJUk8W0ZIkSVJPFtGSJElSTxbRkiRJUk8W0ZIkSVJPFtGSJElSTxbRkiRJUk8W0ZIkSVJPFtGSJElSTxbRkiRJUk8W0ZIkSVJPFtGSJElSTxbRkiRJUk8W0ZIkSVJPFtGSJElSTxbRkiRJUk+diugkO5LcneRQkqsXWf+OJAeT3JHkC0mePbDu8SS3t4/ZYQYvSTqeOVuSRm/jch2SbACuBS4G5oH9SWar6uBAt28AM1X1cJK3Au8DfrVd96OqumDIcUuSFmHOlqTx6HIk+iLgUFUdrqpHgRuAnYMdqurWqnq4XdwHbB1umJKkjszZkjQGXYroM4F7B5bn27alXAl8fmD5tCRzSfYlee1iGyTZ3faZO3r0aIeQJElLGHnOBvO2JC07nQPIIm21aMfkDcAM8PKB5rOr6kiSc4EvJvlmVX3nSS9WtQfYAzAzM7Poa0uSOhl5zgbztiR1ORI9D5w1sLwVOLKwU5JXAe8CLquqR461V9WR9udh4DbgwlXEK0k6MXO2JI1BlyJ6P7A9yTlJTgF2AU86YzvJhcB1NMn4/oH2TUlObZ9vBl4GDJ7cIkkaLnO2JI3BstM5quqxJFcBNwMbgL1VdVeSa4C5qpoFfh/4KeAzSQC+X1WXAecB1yV5gqZgf++CM8QlSUNkzpak8UjVyTWVbWZmpubm5iYdhiStSJIDVTUz6TjGybwtaVqtJmd7x0JJkiSpJ4toSZIkqSeLaEmSJKkni2hJkiSpJ4toSZIkqSeLaEmSJKkni2hJkiSpJ4toSZIkqSeLaEmSJKkni2hJkiSpJ4toSZIkqSeLaEmSJKkni2hJkiSpJ4toSZIkqSeLaEmSJKkni2hJkiSpp05FdJIdSe5OcijJ1YusPzXJp9r1X0uybWDdO9v2u5O8enihS5IWY86WpNFbtohOsgG4FngNcD5weZLzF3S7Eniwqn4W+ADwe+225wO7gJ8DdgB/2L6eJGkEzNmSNB5djkRfBByqqsNV9ShwA7BzQZ+dwMfb5zcCr0yStv2Gqnqkqr4LHGpfT5I0GuZsSRqDjR36nAncO7A8D7xkqT5V9ViSh4Bntu37Fmx75sJfkGQ3sLtdfCTJnZ2iXzs2Az+YdBBj5pjXh/U45udN+PePPGeDeZv1995eb+MFx7xerDhndymis0hbdezTZVuqag+wByDJXFXNdIhrzXDM64NjXh+SzE06hEXahpqzwby93sa83sYLjnm9WE3O7jKdYx44a2B5K3BkqT5JNgJPBx7ouK0kaXjM2ZI0Bl2K6P3A9iTnJDmF5qST2QV9ZoEr2uevA75YVdW272rPBD8H2A58fTihS5IWYc6WpDFYdjpHO1/uKuBmYAOwt6ruSnINMFdVs8BHgU8mOURzNGNXu+1dST4NHAQeA95WVY8v8yv3rHw4U8sxrw+OeX2Y6JgnkLPB/bwerLfxgmNeL1Y85jQHHyRJkiR15R0LJUmSpJ4soiVJkqSeJlZEr+a2tNOqw5jfkeRgkjuSfCHJsycR5zAtN+aBfq9LUkmm/tI6Xcac5FfafX1Xkj8Zd4zD1uG9fXaSW5N8o31/XzqJOIclyd4k9y91beQ0PtT+e9yR5EXjjnHYzNnm7AX9zNlTzJx93PqV5eyqGvuD5mSX7wDnAqcAfwmcv6DPbwAfaZ/vAj41iVjHPOZfBJ7aPn/rehhz2+9pwJdobvIwM+m4x7CftwPfADa1yz8z6bjHMOY9wFvb5+cD35t03Ksc8y8ALwLuXGL9pcDnaa67/FLga5OOeQz72Jy9Dsbc9jNnnwSxj3jM5uwOrzupI9GruS3ttFp2zFV1a1U93C7uo7lG6zTrsp8B3gO8D/j7cQY3Il3G/OvAtVX1IEBV3T/mGIety5gL+On2+dOZ8msPV9WXaK5qsZSdwCeqsQ84I8mzxhPdSJizzdmDzNnTzZx9vBXl7EkV0YvdlnbhrWWfdFta4NhtaadVlzEPupLmf0XTbNkxJ7kQOKuqPjfOwEaoy35+LvDcJF9Osi/JjrFFNxpdxvw7wBuSzAM3AW8fT2gT0/fv/WRnzjZnA+Zsc/aataKc3eW236OwmtvSTqvO40nyBmAGePlIIxq9E445yU8AHwDeNK6AxqDLft5I8/XgK2iOXP33JC+oqh+OOLZR6TLmy4GPVdV/SvLzNNcofkFVPTH68CZiPeav9TjmpqM5e5qZsxvm7OMtm78mdSR6NbelnVadbqeb5FXAu4DLquqRMcU2KsuN+WnAC4DbknyPZh7S7JSfqNL1vf1nVfXjqvoucDdNgp5WXcZ8JfBpgKr6KnAasHks0U3GWrt9tjnbnA3mbHP22rWinD2pIno1t6WdVsuOuf2a7DqaZDztc65gmTFX1UNVtbmqtlXVNpo5hZdV1dxkwh2KLu/tz9KckESSzTRfFR4ea5TD1WXM3wdeCZDkPJqEfHSsUY7XLPDG9ozvlwIPVdV9kw5qFczZ5mxzNubssUY5XivL2RM8U/JS4Ns0Z4i+q227huYPEpod9hngEPB14NxJxTrGMf8F8LfA7e1jdtIxj3rMC/rexpSf6d1xPwd4P82tlb8J7Jp0zGMY8/nAl2nOAr8duGTSMa9yvNcD9wE/pjmCcSXwFuAtA/v42vbf45vr5H1tzjZnT+XDnG3OXmnO9rbfkiRJUk/esVCSJEnqySJakiRJ6skiWpIkSerJIlqSJEnqySJakiRJ6skiWpIkSerJIlqSJEnq6f8C+qGnsgCg/acAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_trace(trace);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NUTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import absolute_import\n",
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "\n",
    "import functools\n",
    "import tensorflow as tf\n",
    "import tensorflow_probability as tfp\n",
    "\n",
    "tfd = tfp.distributions\n",
    "\n",
    "__all__ = [\n",
    "    \"kernel\",\n",
    "]\n",
    "\n",
    "\n",
    "def nuts(target_log_prob_fn,\n",
    "           current_state,\n",
    "           step_size,\n",
    "           seed=None,\n",
    "           current_target_log_prob=None,\n",
    "           current_grads_target_log_prob=None,\n",
    "           name=None):\n",
    "  \"\"\"Simulates a No-U-Turn Sampler (NUTS) trajectory.\n",
    "  Args:\n",
    "    target_log_prob_fn: Python callable which takes an argument like\n",
    "      `*current_state` and returns its (possibly unnormalized) log-density under\n",
    "      the target distribution.\n",
    "    current_state: List of `Tensor`s representing the states to simulate from.\n",
    "    step_size: List of `Tensor`s representing the step sizes for the leapfrog\n",
    "      integrator. Must have same shape as `current_state`.\n",
    "    seed: Integer to seed the random number generator.\n",
    "    current_target_log_prob: Scalar `Tensor` representing the value of\n",
    "      `target_log_prob_fn` at the `current_state`.\n",
    "    current_grads_target_log_prob: List of `Tensor`s representing gradient of\n",
    "      `current_target_log_prob` with respect to `current_state`. Must have same\n",
    "      shape as `current_state`.\n",
    "    name: A name for the operation.\n",
    "  Returns:\n",
    "    next_state: List of `Tensor`s representing the next states of the NUTS\n",
    "      trajectory. Has same shape as `current_state`.\n",
    "    next_target_log_prob: Scalar `Tensor` representing the value of\n",
    "      `target_log_prob_fn` at `next_state`.\n",
    "    next_grads_target_log_prob: List of `Tensor`s representing the gradient of\n",
    "      `next_target_log_prob` with respect to `next_state`.\n",
    "  Raises:\n",
    "    NotImplementedError: If the execution mode is not eager.\n",
    "  \"\"\"\n",
    "  if not tf.executing_eagerly():\n",
    "    raise NotImplementedError(\"`kernel` is only available in Eager mode.\")\n",
    "\n",
    "  with tf.compat.v1.name_scope(\n",
    "      name,\n",
    "      default_name=\"nuts_kernel\",\n",
    "      values=[\n",
    "          current_state, step_size, seed, current_target_log_prob,\n",
    "          current_grads_target_log_prob\n",
    "      ]):\n",
    "    with tf.compat.v1.name_scope(\"initialize\"):\n",
    "      current_state = [tf.convert_to_tensor(value=s) for s in current_state]\n",
    "      step_size = [tf.convert_to_tensor(value=s) for s in step_size]\n",
    "      value_and_gradients_fn = lambda *args: tfp.math.value_and_gradient(  # pylint: disable=g-long-lambda\n",
    "          target_log_prob_fn, args)\n",
    "      value_and_gradients_fn = _embed_no_none_gradient_check(\n",
    "          value_and_gradients_fn)\n",
    "      if (current_target_log_prob is None or\n",
    "          current_grads_target_log_prob is None):\n",
    "        (current_target_log_prob,\n",
    "         current_grads_target_log_prob) = value_and_gradients_fn(*current_state)\n",
    "\n",
    "      seed_stream = tfd.SeedStream(seed, \"nuts_kernel\")\n",
    "      current_momentum = []\n",
    "      for state_tensor in current_state:\n",
    "        momentum_tensor = tf.random.normal(\n",
    "            shape=tf.shape(input=state_tensor),\n",
    "            dtype=state_tensor.dtype,\n",
    "            seed=seed_stream())\n",
    "        current_momentum.append(momentum_tensor)\n",
    "\n",
    "      # Draw a slice variable u ~ Uniform(0, p(initial state, initial\n",
    "      # momentum)) and compute log u. For numerical stability, we perform this\n",
    "      # in log space where log u = log (u' * p(...)) = log u' + log\n",
    "      # p(...) and u' ~ Uniform(0, 1).\n",
    "      log_slice_sample = tf.math.log(tf.random.uniform([], seed=seed_stream()))\n",
    "      log_slice_sample += _log_joint(current_target_log_prob,\n",
    "                                     current_momentum)\n",
    "\n",
    "      # Initialize loop variables. It comprises a collection of information\n",
    "      # about a \"reverse\" state, a collection of information about a \"forward\"\n",
    "      # state, a collection of information about the next state,\n",
    "      # the trajectory's tree depth, the number of candidate states, and\n",
    "      # whether to continue the trajectory.\n",
    "      reverse_state = current_state\n",
    "      reverse_target_log_prob = current_target_log_prob\n",
    "      reverse_grads_target_log_prob = current_grads_target_log_prob\n",
    "      reverse_momentum = current_momentum\n",
    "      forward_state = current_state\n",
    "      forward_target_log_prob = current_target_log_prob\n",
    "      forward_grads_target_log_prob = current_grads_target_log_prob\n",
    "      forward_momentum = current_momentum\n",
    "      next_state = current_state\n",
    "      next_target_log_prob = current_target_log_prob\n",
    "      next_grads_target_log_prob = current_grads_target_log_prob\n",
    "      depth = 0\n",
    "      num_states = 1\n",
    "      continue_trajectory = True\n",
    "\n",
    "    while continue_trajectory:\n",
    "      # Grow the No-U-Turn Sampler trajectory by choosing a random direction and\n",
    "      # simulating Hamiltonian dynamics in that direction. This extends either\n",
    "      # the forward or reverse state.\n",
    "      direction = tfp.math.random_rademacher([], seed=seed_stream())\n",
    "      if direction < 0:\n",
    "        [\n",
    "            reverse_state,\n",
    "            reverse_target_log_prob,\n",
    "            reverse_grads_target_log_prob,\n",
    "            reverse_momentum,\n",
    "            _,\n",
    "            _,\n",
    "            _,\n",
    "            _,\n",
    "            next_state_in_subtree,\n",
    "            next_target_log_prob_in_subtree,\n",
    "            next_grads_target_log_prob_in_subtree,\n",
    "            num_states_in_subtree,\n",
    "            continue_trajectory,\n",
    "        ] = _build_tree(\n",
    "            value_and_gradients_fn=value_and_gradients_fn,\n",
    "            current_state=reverse_state,\n",
    "            current_target_log_prob=reverse_target_log_prob,\n",
    "            current_grads_target_log_prob=reverse_grads_target_log_prob,\n",
    "            current_momentum=reverse_momentum,\n",
    "            direction=direction,\n",
    "            depth=depth,\n",
    "            step_size=step_size,\n",
    "            log_slice_sample=log_slice_sample,\n",
    "            seed=seed_stream())\n",
    "      else:\n",
    "        [\n",
    "            _,\n",
    "            _,\n",
    "            _,\n",
    "            _,\n",
    "            forward_state,\n",
    "            forward_target_log_prob,\n",
    "            forward_grads_target_log_prob,\n",
    "            forward_momentum,\n",
    "            next_state_in_subtree,\n",
    "            next_target_log_prob_in_subtree,\n",
    "            next_grads_target_log_prob_in_subtree,\n",
    "            num_states_in_subtree,\n",
    "            continue_trajectory,\n",
    "        ] = _build_tree(\n",
    "            value_and_gradients_fn=value_and_gradients_fn,\n",
    "            current_state=forward_state,\n",
    "            current_target_log_prob=forward_target_log_prob,\n",
    "            current_grads_target_log_prob=forward_grads_target_log_prob,\n",
    "            current_momentum=forward_momentum,\n",
    "            direction=direction,\n",
    "            depth=depth,\n",
    "            step_size=step_size,\n",
    "            log_slice_sample=log_slice_sample,\n",
    "            seed=seed_stream())\n",
    "\n",
    "      if continue_trajectory:\n",
    "        # If the built tree did not terminate, accept the tree's next state\n",
    "        # with a certain probability.\n",
    "        accept_state_in_subtree = _random_bernoulli(\n",
    "            [],\n",
    "            probs=tf.minimum(1., num_states_in_subtree / num_states),\n",
    "            dtype=tf.bool,\n",
    "            seed=seed_stream())\n",
    "        if accept_state_in_subtree:\n",
    "          next_state = next_state_in_subtree\n",
    "          next_target_log_prob = next_target_log_prob_in_subtree\n",
    "          next_grads_target_log_prob = next_grads_target_log_prob_in_subtree\n",
    "\n",
    "      # Continue the NUTS trajectory if the tree-building did not terminate, and\n",
    "      # if the reverse-most and forward-most states do not exhibit a U-turn.\n",
    "      has_no_u_turn = tf.logical_and(\n",
    "          _has_no_u_turn(forward_state, reverse_state, forward_momentum),\n",
    "          _has_no_u_turn(forward_state, reverse_state, reverse_momentum))\n",
    "      continue_trajectory = continue_trajectory and has_no_u_turn\n",
    "      num_states += num_states_in_subtree\n",
    "      depth += 1\n",
    "\n",
    "    return next_state, next_target_log_prob, next_grads_target_log_prob\n",
    "\n",
    "\n",
    "def _build_tree(value_and_gradients_fn,\n",
    "                current_state,\n",
    "                current_target_log_prob,\n",
    "                current_grads_target_log_prob,\n",
    "                current_momentum,\n",
    "                direction,\n",
    "                depth,\n",
    "                step_size,\n",
    "                log_slice_sample,\n",
    "                max_simulation_error=1000.,\n",
    "                seed=None):\n",
    "  \"\"\"Builds a tree at a given tree depth and at a given state.\n",
    "  The `current` state is immediately adjacent to, but outside of,\n",
    "  the subtrajectory spanned by the returned `forward` and `reverse` states.\n",
    "  Args:\n",
    "    value_and_gradients_fn: Python callable which takes an argument like\n",
    "      `*current_state` and returns a tuple of its (possibly unnormalized)\n",
    "      log-density under the target distribution and its gradient with respect to\n",
    "      each state.\n",
    "    current_state: List of `Tensor`s representing the current states of the\n",
    "      NUTS trajectory.\n",
    "    current_target_log_prob: Scalar `Tensor` representing the value of\n",
    "      `target_log_prob_fn` at the `current_state`.\n",
    "    current_grads_target_log_prob: List of `Tensor`s representing gradient of\n",
    "      `current_target_log_prob` with respect to `current_state`. Must have same\n",
    "      shape as `current_state`.\n",
    "    current_momentum: List of `Tensor`s representing the momentums of\n",
    "      `current_state`. Must have same shape as `current_state`.\n",
    "    direction: int that is either -1 or 1. It determines whether to perform\n",
    "      leapfrog integration backwards (reverse) or forward in time respectively.\n",
    "    depth: non-negative int that indicates how deep of a tree to build.\n",
    "      Each call to `_build_tree` takes `2**depth` leapfrog steps.\n",
    "    step_size: List of `Tensor`s representing the step sizes for the leapfrog\n",
    "      integrator. Must have same shape as `current_state`.\n",
    "    log_slice_sample: The log of an auxiliary slice variable. It is used\n",
    "      together with `max_simulation_error` to avoid simulating trajectories with\n",
    "      too much numerical error.\n",
    "    max_simulation_error: Maximum simulation error to tolerate before\n",
    "      terminating the trajectory. Simulation error is the\n",
    "      `log_slice_sample` minus the log-joint probability at the simulated state.\n",
    "    seed: Integer to seed the random number generator.\n",
    "  Returns:\n",
    "    reverse_state: List of `Tensor`s representing the \"reverse\" states of the\n",
    "      NUTS trajectory. Has same shape as `current_state`.\n",
    "    reverse_target_log_prob: Scalar `Tensor` representing the value of\n",
    "      `target_log_prob_fn` at the `reverse_state`.\n",
    "    reverse_grads_target_log_prob: List of `Tensor`s representing gradient of\n",
    "      `reverse_target_log_prob` with respect to `reverse_state`. Has same shape\n",
    "      as `reverse_state`.\n",
    "    reverse_momentum: List of `Tensor`s representing the momentums of\n",
    "      `reverse_state`. Has same shape as `reverse_state`.\n",
    "    forward_state: List of `Tensor`s representing the \"forward\" states of the\n",
    "      NUTS trajectory. Has same shape as `current_state`.\n",
    "    forward_target_log_prob: Scalar `Tensor` representing the value of\n",
    "      `target_log_prob_fn` at the `forward_state`.\n",
    "    forward_grads_target_log_prob: List of `Tensor`s representing gradient of\n",
    "      `forward_target_log_prob` with respect to `forward_state`. Has same shape\n",
    "      as `forward_state`.\n",
    "    forward_momentum: List of `Tensor`s representing the momentums of\n",
    "      `forward_state`. Has same shape as `forward_state`.\n",
    "    next_state: List of `Tensor`s representing the next states of the NUTS\n",
    "      trajectory. Has same shape as `current_state`.\n",
    "    next_target_log_prob: Scalar `Tensor` representing the value of\n",
    "      `target_log_prob_fn` at `next_state`.\n",
    "    next_grads_target_log_prob: List of `Tensor`s representing the gradient of\n",
    "      `next_target_log_prob` with respect to `next_state`.\n",
    "    num_states: Number of acceptable candidate states in the subtree. A state is\n",
    "      acceptable if it is \"in the slice\", that is, if its log-joint probability\n",
    "      with its momentum is greater than `log_slice_sample`.\n",
    "    continue_trajectory: bool determining whether to continue the simulation\n",
    "      trajectory. The trajectory is continued if no U-turns are encountered\n",
    "      within the built subtree, and if the log-probability accumulation due to\n",
    "      integration error does not exceed `max_simulation_error`.\n",
    "  \"\"\"\n",
    "  if depth == 0:  # base case\n",
    "    # Take a leapfrog step. Terminate the tree-building if the simulation\n",
    "    # error from the leapfrog integrator is too large. States discovered by\n",
    "    # continuing the simulation are likely to have very low probability.\n",
    "    [\n",
    "        next_state,\n",
    "        next_target_log_prob,\n",
    "        next_grads_target_log_prob,\n",
    "        next_momentum,\n",
    "    ] = _leapfrog(\n",
    "        value_and_gradients_fn=value_and_gradients_fn,\n",
    "        current_state=current_state,\n",
    "        current_grads_target_log_prob=current_grads_target_log_prob,\n",
    "        current_momentum=current_momentum,\n",
    "        step_size=direction * step_size)\n",
    "    next_log_joint = _log_joint(next_target_log_prob, next_momentum)\n",
    "    num_states = tf.cast(next_log_joint > log_slice_sample, dtype=tf.int32)\n",
    "    continue_trajectory = (next_log_joint >\n",
    "                           log_slice_sample - max_simulation_error)\n",
    "    return [\n",
    "        next_state,\n",
    "        next_target_log_prob,\n",
    "        next_grads_target_log_prob,\n",
    "        next_momentum,\n",
    "        next_state,\n",
    "        next_target_log_prob,\n",
    "        next_grads_target_log_prob,\n",
    "        next_momentum,\n",
    "        next_state,\n",
    "        next_target_log_prob,\n",
    "        next_grads_target_log_prob,\n",
    "        num_states,\n",
    "        continue_trajectory,\n",
    "    ]\n",
    "\n",
    "  # Build a tree at the current state.\n",
    "  seed_stream = tfd.SeedStream(seed, \"build_tree\")\n",
    "  [\n",
    "      reverse_state,\n",
    "      reverse_target_log_prob,\n",
    "      reverse_grads_target_log_prob,\n",
    "      reverse_momentum,\n",
    "      forward_state,\n",
    "      forward_target_log_prob,\n",
    "      forward_grads_target_log_prob,\n",
    "      forward_momentum,\n",
    "      next_state,\n",
    "      next_target_log_prob,\n",
    "      next_grads_target_log_prob,\n",
    "      num_states,\n",
    "      continue_trajectory,\n",
    "  ] = _build_tree(value_and_gradients_fn=value_and_gradients_fn,\n",
    "                  current_state=current_state,\n",
    "                  current_target_log_prob=current_target_log_prob,\n",
    "                  current_grads_target_log_prob=current_grads_target_log_prob,\n",
    "                  current_momentum=current_momentum,\n",
    "                  direction=direction,\n",
    "                  depth=depth - 1,\n",
    "                  step_size=step_size,\n",
    "                  log_slice_sample=log_slice_sample,\n",
    "                  seed=seed_stream())\n",
    "  if continue_trajectory:\n",
    "    # If the just-built subtree did not terminate, build a second subtree at\n",
    "    # the forward or reverse state, as appropriate.\n",
    "    if direction < 0:\n",
    "      [\n",
    "          reverse_state,\n",
    "          reverse_target_log_prob,\n",
    "          reverse_grads_target_log_prob,\n",
    "          reverse_momentum,\n",
    "          _,\n",
    "          _,\n",
    "          _,\n",
    "          _,\n",
    "          far_state,\n",
    "          far_target_log_prob,\n",
    "          far_grads_target_log_prob,\n",
    "          far_num_states,\n",
    "          far_continue_trajectory,\n",
    "      ] = _build_tree(\n",
    "          value_and_gradients_fn=value_and_gradients_fn,\n",
    "          current_state=reverse_state,\n",
    "          current_target_log_prob=reverse_target_log_prob,\n",
    "          current_grads_target_log_prob=reverse_grads_target_log_prob,\n",
    "          current_momentum=reverse_momentum,\n",
    "          direction=direction,\n",
    "          depth=depth - 1,\n",
    "          step_size=step_size,\n",
    "          log_slice_sample=log_slice_sample,\n",
    "          seed=seed_stream())\n",
    "    else:\n",
    "      [\n",
    "          _,\n",
    "          _,\n",
    "          _,\n",
    "          _,\n",
    "          forward_state,\n",
    "          forward_target_log_prob,\n",
    "          forward_grads_target_log_prob,\n",
    "          forward_momentum,\n",
    "          far_state,\n",
    "          far_target_log_prob,\n",
    "          far_grads_target_log_prob,\n",
    "          far_num_states,\n",
    "          far_continue_trajectory,\n",
    "      ] = _build_tree(\n",
    "          value_and_gradients_fn=value_and_gradients_fn,\n",
    "          current_state=forward_state,\n",
    "          current_target_log_prob=forward_target_log_prob,\n",
    "          current_grads_target_log_prob=forward_grads_target_log_prob,\n",
    "          current_momentum=forward_momentum,\n",
    "          direction=direction,\n",
    "          depth=depth - 1,\n",
    "          step_size=step_size,\n",
    "          log_slice_sample=log_slice_sample,\n",
    "          seed=seed_stream())\n",
    "\n",
    "    # Propose either `next_state` (which came from the first subtree and so is\n",
    "    # nearby) or the new forward/reverse state (which came from the second\n",
    "    # subtree and so is far away).\n",
    "    num_states += far_num_states\n",
    "    accept_far_state = _random_bernoulli(\n",
    "        [],\n",
    "        probs=far_num_states / num_states,\n",
    "        dtype=tf.bool,\n",
    "        seed=seed_stream())\n",
    "    if accept_far_state:\n",
    "      next_state = far_state\n",
    "      next_target_log_prob = far_target_log_prob\n",
    "      next_grads_target_log_prob = far_grads_target_log_prob\n",
    "\n",
    "    # Continue the NUTS trajectory if the far subtree did not terminate either,\n",
    "    # and if the reverse-most and forward-most states do not exhibit a U-turn.\n",
    "    has_no_u_turn = tf.logical_and(\n",
    "        _has_no_u_turn(forward_state, reverse_state, forward_momentum),\n",
    "        _has_no_u_turn(forward_state, reverse_state, reverse_momentum))\n",
    "    continue_trajectory = far_continue_trajectory and has_no_u_turn\n",
    "\n",
    "  return [\n",
    "      reverse_state,\n",
    "      reverse_target_log_prob,\n",
    "      reverse_grads_target_log_prob,\n",
    "      reverse_momentum,\n",
    "      forward_state,\n",
    "      forward_target_log_prob,\n",
    "      forward_grads_target_log_prob,\n",
    "      forward_momentum,\n",
    "      next_state,\n",
    "      next_target_log_prob,\n",
    "      next_grads_target_log_prob,\n",
    "      num_states,\n",
    "      continue_trajectory,\n",
    "  ]\n",
    "\n",
    "\n",
    "def _embed_no_none_gradient_check(value_and_gradients_fn):\n",
    "  \"\"\"Wraps value and gradients function to assist with None gradients.\"\"\"\n",
    "  @functools.wraps(value_and_gradients_fn)\n",
    "  def func_wrapped(*args, **kwargs):\n",
    "    \"\"\"Wrapped function which checks for None gradients.\"\"\"\n",
    "    value, grads = value_and_gradients_fn(*args, **kwargs)\n",
    "    if any(grad is None for grad in grads):\n",
    "      raise ValueError(\"Gradient is None for a state.\")\n",
    "    return value, grads\n",
    "  return func_wrapped\n",
    "\n",
    "\n",
    "def _has_no_u_turn(state_one, state_two, momentum):\n",
    "  \"\"\"If two given states and momentum do not exhibit a U-turn pattern.\"\"\"\n",
    "  dot_product = sum([\n",
    "      tf.reduce_sum(input_tensor=(s1 - s2) * m)\n",
    "      for s1, s2, m in zip(state_one, state_two, momentum)\n",
    "  ])\n",
    "  return dot_product > 0\n",
    "\n",
    "\n",
    "def _leapfrog(value_and_gradients_fn,\n",
    "              current_state,\n",
    "              current_grads_target_log_prob,\n",
    "              current_momentum,\n",
    "              step_size):\n",
    "  \"\"\"Runs one step of leapfrog integration.\"\"\"\n",
    "  mid_momentum = [\n",
    "      m + 0.5 * step * g for m, step, g in\n",
    "      zip(current_momentum, step_size, current_grads_target_log_prob)]\n",
    "  next_state = [\n",
    "      s + step * m for s, step, m in\n",
    "      zip(current_state, step_size, mid_momentum)]\n",
    "  next_target_log_prob, next_grads_target_log_prob = value_and_gradients_fn(\n",
    "      *next_state)\n",
    "  next_momentum = [\n",
    "      m + 0.5 * step * g for m, step, g in\n",
    "      zip(mid_momentum, step_size, next_grads_target_log_prob)]\n",
    "  return [\n",
    "      next_state,\n",
    "      next_target_log_prob,\n",
    "      next_grads_target_log_prob,\n",
    "      next_momentum,\n",
    "  ]\n",
    "\n",
    "\n",
    "def _log_joint(current_target_log_prob, current_momentum):\n",
    "  \"\"\"Log-joint probability given a state's log-probability and momentum.\"\"\"\n",
    "  momentum_log_prob = -sum(\n",
    "      [tf.reduce_sum(input_tensor=0.5 * (m**2.)) for m in current_momentum])\n",
    "  return current_target_log_prob + momentum_log_prob\n",
    "\n",
    "\n",
    "def _random_bernoulli(shape, probs, dtype=tf.int32, seed=None, name=None):\n",
    "  \"\"\"Returns samples from a Bernoulli distribution.\"\"\"\n",
    "  with tf.compat.v1.name_scope(name, \"random_bernoulli\", [shape, probs]):\n",
    "    probs = tf.convert_to_tensor(value=probs)\n",
    "    random_uniform = tf.random.uniform(shape, dtype=probs.dtype, seed=seed)\n",
    "    return tf.cast(tf.less(random_uniform, probs), dtype)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_nuts(model, nsteps=1000, burnin=500, step_size=.05, profile=False, compile=False):\n",
    "    # Since HMC operates over unconstrained space, we need to transform the\n",
    "    # samples so they live in real-space.\n",
    "    random_variables = model._forward_context.vars\n",
    "    unnormalized_posterior_log_prob = model.make_log_prob_function()\n",
    "    inits = [0.] * len(random_variables)\n",
    "    if compile:\n",
    "        unnormalized_posterior_log_prob = tf.function(\n",
    "            unnormalized_posterior_log_prob)\n",
    "\n",
    "    posterior_samples = []\n",
    "    target_log_prob = None\n",
    "    grads_target_log_prob = None\n",
    "    \n",
    "    if profile:\n",
    "        sampler = profiler(nuts)\n",
    "    else:\n",
    "        sampler = nuts\n",
    "        \n",
    "    step_size = tf.stack([step_size] * len(random_variables))\n",
    "    \n",
    "    for step in range(nsteps):\n",
    "        if step % 10 == 0:\n",
    "            print(\"Step\", step)\n",
    "        [\n",
    "            inits,\n",
    "            target_log_prob,\n",
    "            grads_target_log_prob,\n",
    "        ] = sampler(target_log_prob_fn=unnormalized_posterior_log_prob,\n",
    "                   current_state=inits,\n",
    "                   step_size=step_size,\n",
    "                   seed=step,\n",
    "                   current_target_log_prob=target_log_prob,\n",
    "                   current_grads_target_log_prob=grads_target_log_prob)\n",
    "        posterior_samples.append(inits)\n",
    "   \n",
    "    trace = {rv.name: arr for rv, arr in zip(random_variables, posterior_samples)}\n",
    "    \n",
    "    return az.dict_to_dataset(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 0\n",
      "Step 10\n",
      "Step 20\n",
      "Step 30\n",
      "Step 40\n",
      "Step 50\n",
      "Step 60\n",
      "Step 70\n",
      "Step 80\n",
      "Step 90\n",
      "CPU times: user 16.9 s, sys: 82.4 ms, total: 17 s\n",
      "Wall time: 17.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "trace = sample_nuts(model, 100, profile=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with defun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 0\n",
      "Step 10\n",
      "Step 20\n",
      "Step 30\n",
      "Step 40\n",
      "Step 50\n",
      "Step 60\n",
      "Step 70\n",
      "Step 80\n",
      "Step 90\n",
      "CPU times: user 6.65 s, sys: 214 ms, total: 6.87 s\n",
      "Wall time: 6.59 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "trace = sample_nuts(model, 100, compile=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "histogram(): `bins` should be a positive integer",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-53-7c3f6fe9d699>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0maz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot_trace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/arviz/plots/traceplot.py\u001b[0m in \u001b[0;36mplot_trace\u001b[0;34m(data, var_names, coords, divergences, figsize, textsize, lines, combined, kde_kwargs, hist_kwargs, trace_kwargs)\u001b[0m\n\u001b[1;32m    148\u001b[0m                 \u001b[0m_histplot_op\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0midx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhist_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    149\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 150\u001b[0;31m                 \u001b[0mplot_kde\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtextsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mxt_labelsize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0midx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkde_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    151\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    152\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkind\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"i\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/arviz/plots/kdeplot.py\u001b[0m in \u001b[0;36mplot_kde\u001b[0;34m(values, values2, cumulative, rug, label, bw, quantiles, rotated, contour, fill_last, textsize, plot_kwargs, fill_kwargs, rug_kwargs, contour_kwargs, ax, legend)\u001b[0m\n\u001b[1;32m    168\u001b[0m         \u001b[0mrug_kwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msetdefault\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"markersize\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mmarkersize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    169\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 170\u001b[0;31m         \u001b[0mdensity\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlower\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupper\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_fast_kde\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcumulative\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    171\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    172\u001b[0m         \u001b[0mrug_space\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdensity\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mrug_kwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"space\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/arviz/plots/kdeplot.py\u001b[0m in \u001b[0;36m_fast_kde\u001b[0;34m(x, cumulative, bw, xmin, xmax)\u001b[0m\n\u001b[1;32m    283\u001b[0m     \u001b[0mn_bins\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen_x\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mstd_x\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_points\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    284\u001b[0m     \u001b[0md_x\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mxmax\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mxmin\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mn_bins\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 285\u001b[0;31m     \u001b[0mgrid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_histogram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_bins\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrange_hist\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxmin\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxmax\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    286\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    287\u001b[0m     \u001b[0mscotts_factor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen_x\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m0.2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: histogram(): `bins` should be a positive integer"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEoCAYAAAAt0dJ4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X1snPd16Pnv4etIJIe2RIpULCtyruU4rlFsUq6TbnBbF04Kx1hY/ySpXQRJuk60zb3uBZrcoL5IkQYuFnCT7Qbtrm9TNTHyctHYTrCbCL3u9W7zgnSLOLCM3GZjF95VXddmbZGy7JB6G77N2T9mrND0UBxSM5yH4vcDEJ5nnp9nDn4Y8ujM7zy/JzITSZIkSVLndXU6AEmSJElSjQWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFcSaBVpEPBAR0xHx01XOR0T8aUQcj4ifRMTbWh+mJEmSJF3+mllB+zJw60XOvwc4WP85DPzZpYclSZIkSdvPmgVaZv4AePkiQw4BX82ax4ArImJvqwKUJEmSpO2iFdegXQU8v+x4sv6cJEmSJGkdelrwGtHguWw4MOIwtTZIBgYGfun6669vwdtLkraqJ5544qXMHO10HK02MjKSBw4c6HQYkqQOuZT81ooCbRK4etnxPuCFRgMz8whwBGBiYiKPHTvWgreXJG1VEfHPnY6hHQ4cOIA5TpK2r0vJb61ocTwKfLC+m+M7gJnMfLEFrytJkiRJ28qaK2gR8XXgZmAkIiaBPwB6ATLzC8AjwG3AceAc8FvtClaSJEmSLmdrFmiZeeca5xP4ty2LSJIkSZI2UWVhienZOaZOV3jjrp3sKZc6FksrrkGTJEmSpC1l5vwCU7MVpmYrvHx2HoBSbzfjHSzOwAJNkiRJ0jawVE1OnZnjxGyFEzMVzi8sAXDFzj6uHy8zXi4xvLO3w1FaoEmSJEm6TFUWluqrZHNMn66wVE26u4I9QyXeXO5nrFyi1Nvd6TBfwwJNkiRJ0mVj5txCbZVstsLPztVaF3f0drN/107GyyV2D/bT3dXoVs7FYIEmSZIkactaqiYnT88xVS/KKvXWxV0Dfbxlb5mxconhHZ1vXWyWBZokSZKkLaWysMSJmdoGHyfPzLFUTXq6uthT7mdsqMSecn/hWhebZYEmSZIkqdAyk5nzCxc2+Jg5vwDAzr4e3rhrgLHhfkYG+ukqcOtisyzQJEmSJBXO4lKVk2fmmJqd48RMhbnFJSKCXTv7uGFvmbHhEuXS1mldbJYFmiRJkqRCOD+/dGGV7KUzc1Qz6e3uYs9QbcfFPeV++nu2ZutisyzQJEmSJHVEZvKzZbsuztZbFwf6erhmZICxcondA32XRetisyzQJEnbVkTcCvwJ0A18MTPvW3G+H/gq8EvAKeA3MvPZzY5Tki4nC0vVC7suTs3OXWhd3D3Qxy+8YZixcj9Dl2HrYrMs0CRJ21JEdAP3A+8GJoHHI+JoZj61bNhdwCuZeW1E3AH8EfAbmx+tJG1t5+YXOTFTWyU7dWaeaiZ93fVdF8sl9gyV6Ovp6nSYhWCBJknarm4CjmfmMwAR8SBwCFheoB0CPlN//E3gf4uIyMzczEAlaavJTF45t3BhK/zZSq11cahUa10cHy6xa+f2al1slgWaJGm7ugp4ftnxJPD21cZk5mJEzAC7gZdWvlhEHAYOA+zfv78d8UpSoS0sVZmuty5Oz1aYW6wSEYwM9HHjVcOMlUsM9lt+rKWpGWqiR38/8BXgivqYezLzkRbHKklSKzX62nblylgzY2pPZh4BjgBMTEy4wiZpWzg7t8iJ2QpTMxVOnV3eulhifLjE6GC/rYvrtGaB1mSP/u8DD2fmn0XEDcAjwIE2xCtJUqtMAlcvO94HvLDKmMmI6AGGgZc3JzxJKp7M5OWz87WibLbC6coiAOVSL28aHWC8XGLXQB8Rti5uVDMraM306CdQrj8e5vUJTpKkonkcOBgR1wD/AtwB/OaKMUeBDwE/BN4LfNfrzyRtN/OLVaZP13ZcnJ6tML9UpSuC3YN9HNhd2wp/wNbFlmlmJpvp0f8M8H9GxO8AA8C7WhKdJEltUr+m7G7gUWrt+Q9k5pMRcS9wLDOPAl8CvhYRx6mtnN3RuYglafOcmVu8sMHHqbPzZCb9PV2MDZcYL5cYHeqnt9vWxXZopkBrpv/+TuDLmfnHEfHL1JLZjZlZfc0LeQG1JKlA6tdLP7LiuU8ve1wB3rfZcUnSZqtWk5fPzV8oys7M1VsXd/Ry7egg48MlrtzZa+viJmimQGumR/8u4FaAzPxhRJSAEWB6+SAvoJYkSZKKYX6xWttx8XSF6dm5C62LI4P9vGl0kLFyPzv7bF3cbM3MeDM9+s8BtwBfjoi3ACXgZCsDlSRJknRpTlcWmJqtcGJmjpfPvdq62M348M93XeyxdbGj1izQmuzR/wTwFxHxu9TaHz/sRdSSJElSZ1Wryamz8/WirMLZ+Vrr4vCOXg7uGWS8XOIKWxcLpak1yyZ69J8C3tna0CRJkiSt19ziEtOz9RtGn55jod66ODrUz7V7Bhkrl9jR193pMLUKm0olSZKkLW62ssDUTIUTsxVePjsPQKm3mzdcsYPxcomRwT5bF7cICzRJkiRpi6lWk5fOzjE1M8eJ2Qrn6q2LV+zs483jQ4yXSwzvsHVxK7JAkyRJkraAykK9dbG+6+JitUp3VzA62M91Y7XWxVKvrYtbnQWaJEmSVFAz5xcubPDxyrmfty7uu3IH48MlRgb76e5ylexyYoEmSZIkFcRSNTl1pta2eGKmwvmFJQCu3NnH9ePlWuvizt4OR6l2skCTJEmSOqiysMTUbIWp2TmmT1dYqibdXcGeoRLXl0vsKffburiNWKBJkiRJm2zm3EJtlWy2ws/qrYs7ervZv2sn4+USu21d3LYs0CRJkqQ2W6omJ0/X7k12YrZCpd66uGugj7fsLTNW33VRskCTJEmS2qCysMSJmQpTsxVOnpljqZr0dHWxp9zPeLnE6JCti3o9CzRJkiSpBTKTmfMLFzb4mDm/AMDOvh7euGuAseF+Rgb66bJ1URdhgSZJkiRt0OJSlZNn5pianePETIW5xSUigl07+7hhb5mx4RLlkq2Lap4FmiRJkrQO5+eXLqySvXRmjmomvd1d7BnqZ6y+62J/j62L2hgLNEnSthMRu4CHgAPAs8D7M/OVBuP+C/AO4P/OzP9+M2OUVByZySvnFi5s8DFbb10c6OvhmpEBxsoldg/02bqolrBAkyRtR/cA38nM+yLinvrx7zUY9zlgJ/A/bmZwkjpvYal6YdfFqdm5C62Luwf6+IU3DDNW7mfI1kW1QVMFWkTcCvwJ0A18MTPvazDm/cBngAT+PjN/s4VxSpLUSoeAm+uPvwJ8nwYFWmZ+JyJuXvm8pMvTuflFTszUVslOnZmnmklfd23XxbFyiT1DJfp6ujodpi5zaxZoEdEN3A+8G5gEHo+Io5n51LIxB4H/ALwzM1+JiD3tCliSpBYYy8wXATLzxVbkrYg4DBwG2L9//6W+nKRN8Grr4qtb4c9Waq2LQ6Va6+L4cIldO21d1OZqZgXtJuB4Zj4DEBEPUvvm8allYz4K3P9q/35mTrc6UEmS1iMi/gYYb3DqU+14v8w8AhwBmJiYyHa8h6RLt7BUZbreujg9W2FusUpEMDLQx41XDTNWLjHY71VA6pxmPn1XAc8vO54E3r5izHUAEfF31NogP5OZ/6UlEUqStAGZ+a7VzkXEVETsra+e7QX8YlG6jJ2dW+TEbIWpmQovnZ0nL7QulhgfLjE62G/rogqjmQKt0Zruym8Ge4CD1Pr59wF/GxE3ZubPXvNCtn9IkorhKPAh4L76f7/d2XAktVJm8vLZ+VpRNlvhdGURgHKpl381OsB4ucSugT4ibF1U8TRToE0CVy873ge80GDMY5m5APxTRDxNrWB7fPkg2z8kSQVxH/BwRNwFPAe8DyAiJoDfzsyP1I//FrgeGIyISeCuzHy0QzFLuoj5xSrTp2s7Lk7PVphfqtIVwe7BPg7srm2FP2DroraAZj6ljwMHI+Ia4F+AO4CVOzR+C7gT+HJEjFBreXymlYFKktQqmXkKuKXB88eAjyw7/tebGZek9Tkzt3hhg49T9dbF/p4uxoZLjJdLjA7109tt66K2ljULtMxcjIi7gUepXV/2QGY+GRH3Ascy82j93K9HxFPAEvDJevKTJEmSWqJaTV4+N3+hKDszV29d3NHLtaODjA+XuHJnr62L2tKaWufNzEeAR1Y89+lljxP4eP1HkiRJaon5xWptx8XTFaZn5y60Lo4M9vOm0UHGyv3s7LN1UZcPP82SJEkqlNOVBaZmK5yYmePlc6+2LnYzPvzzXRd7bF3UZcoCTZIkSR1VrSanzs7Xi7IKZ+drrYvDO3o5uGeQ8XKJK2xd1DZhgSZJkqRNN7e4xPRs/YbRp+dYqLcujg71c+2eQcbKJXb0dXc6TGnTWaBJkiRpU8xWFpiaqXBitsLLZ+cBKPV284YrdjBeLjEy2GfrorY9CzRJkiS1RbWavHR2jqmZOU7MVjhXb128Ymcfbx4fYrxcYniHrYvSchZokiRJapnKQr11sb7r4mK1SndXMDrYz3VjtdbFUq+ti9JqLNAkSZJ0SWbOL1zY4OOVc7XWxR293ey7cgfjwyVGBvvp7nKVTGqGBZokSZLWZamanDpTa1s8MVPh/MISAFfu7OP68XKtdXFnb4ejlLYmCzRJkiStqbKwxNRshanZOaZPV1iqJj1dXYwO9XN9ucSecr+ti1ILWKBJkiSpoZlzC7VVstkKP1vWurh/107GyyV227ootZwFmiRJkoBa6+LJ07V7k52YrVCpty7uGujjLXvLjNV3XZTUPhZokiRJ21hlYYkTMxWmZiucPDN3oXVxT7mf8XKJ0SFbF6XNZIEmSZK0jWQmM+cXLmzwMXN+AYCBvh7euGuAseF+Rgb66bJ1UeoICzRJkqTL3OJSlZNn5pianePETIW5xSUigl07+7hhb5mx4RLlkq2LUhFYoEmStp2I2AU8BBwAngXen5mvrBjz3wB/BpSBJeB/ysyHNjdSaePOzy9dWCV76cwc1Ux6u7vYM9TPWH3Xxf4eWxelommqQIuIW4E/AbqBL2bmfauMey/wDeC/zcxjLYtSkqTWugf4TmbeFxH31I9/b8WYc8AHM/P/i4g3AE9ExKOZ+bPNDlZqRmbyyrmFCxt8zNZbFwf7e7hmZICxcondA322LkoFt2aBFhHdwP3Au4FJ4PGIOJqZT60YNwT8O+BH7QhUkqQWOgTcXH/8FeD7rCjQMvP/Xfb4hYiYBkYBCzQVxsJS9cKui1OzcxdaF3cP9PELbxhmrNzPkK2L0pbSzAraTcDxzHwGICIepJbYnlox7g+BzwL/vqURSpLUemOZ+SJAZr4YEXsuNjgibgL6gH+8yJjDwGGA/fv3tzBU6bXOzS9yYqa2SnbqzDzVTPq6a7sujpVL7Bkq0dfT1ekwJW1QMwXaVcDzy44ngbcvHxARbwWuzsy/iohVCzSTlyRps0TE3wDjDU59ap2vsxf4GvChzKyuNi4zjwBHACYmJnI97yFdTGby8tl5pmZrK2WzlVrr4lCp1ro4Plxi105bF6XLRTMFWqPf9guJJyK6gM8DH17rhUxekqTNkpnvWu1cRExFxN766tleYHqVcWXgPwO/n5mPtSlU6XUWlqpMn67tuHjydIW5xSoRwchAHzdeNcxYucRgv3u9SZejZn6zJ4Grlx3vA15YdjwE3Ah8PyKg9m3l0Yi43Y1CJEkFdRT4EHBf/b/fXjkgIvqA/wP4amZ+Y3PD03Z0dm6RE7MVpmYqvHR2nrzQulhifLjE6GC/rYvSNtBMgfY4cDAirgH+BbgD+M1XT2bmDDDy6nFEfB/49xZnkqQCuw94OCLuAp4D3gcQERPAb2fmR4D3A78C7I6ID9f/vw9n5n/tQLy6DGUmp87O1zf4qHC6sghAudTLvxodYLxcYtdAH/UvwCVtE2sWaJm5GBF3A49S22b/gcx8MiLuBY5l5tF2BylJUitl5inglgbPHwM+Un/8n4D/tMmh6TI3v1hl+nStIJuenWN+qUpXBLsH+ziwu7YV/oCti9K21tRfgMx8BHhkxXOfXmXszZceliRJ0uXhzFxt18Wp2Qqn6q2L/T1djA2XGC+XGB3qp7fb1kVJNX5FI0mS1ELV6mtbF8/M1VsXd/Ry7egg48MlrtzZa+uipIYs0CRJki7R3OIS0/Vt8KdPz7FQb10cGeznTaODjJX72dnnP7skrc2/FJIkSRtwurLA1GyFEzNzvHzu1dbFbvYO/3zXxR5bFyWtkwWaJElSE6rV5KWzc0zP1u5Pdna+1ro4vKOXg3sGGS+XuMLWRUmXyAJNkiRpFa+2Lp6YrXByWevi6FA/1+4ZZKxcYkdfd6fDlHQZsUCTJElaZraywNRMhROzFV4+Ow9AqbebN1yxg/FyiZHBPlsXJbWNBZokSdrWqtXkpTNzTNVXys7VWxev2NnHm8eHGC+XGN5h66KkzWGBJkmStp3KwmtbFxerVbq7gtHBfq4bq7UulnptXZS0+SzQJEnStjBz/tVdFyu8cq7Wurijt5t9V+5gfLjEyGA/3V2ukknqLAs0SZJ0WVq60LpYK8rOLywBcOXOPq4fL9daF3f2djhKSXotCzRJknTZqCwsXSjITp6ZY6ma9HR1MTrUz/XlEnvK/bYuSio0CzRJkrSlzZxb4MRsbdfFny1rXdy/ayfj5RK7bV2UtIVYoEmSpC1lqZqcPF1vXZytUKm3Lu4a6OMte8uM1XddlKStyAJNkiQVXmVhiRP1e5O9tKx1cU+5n/FyidEhWxclXR6aKtAi4lbgT4Bu4IuZed+K8x8HPgIsAieB/yEz/7nFsUqSpG0iM5k5X29dnKkwc34BgIG+Ht64a4Cx4X5GBvrpsnVR0mVmzQItIrqB+4F3A5PA4xFxNDOfWjbsx8BEZp6LiI8BnwV+ox0BS5LUChGxC3gIOAA8C7w/M19ZMeaNwP9O7QvKXuB/zcwvbG6k28fiUpWTr94weqbC3OISEcGunX3csLfM2HCJcsnWRUmXt2ZW0G4CjmfmMwAR8SBwCLhQoGXm95aNfwz4QCuDlCSpDe4BvpOZ90XEPfXj31sx5kXgv8vMuYgYBH5a/5Lyhc0O9nJ1fn7pwirZS2fmqGbS293FnqF+xuq7Lvb32LooaftopkC7Cnh+2fEk8PaLjL8L+OtGJyLiMHAYYP/+/U2GKElSWxwCbq4//grwfVYUaJk5v+ywH+jajMAuZ5nJK+cWLmzwMVtvXRzs7+GakQHGyiV2D/TZuihp22qmQGv0FzIbDoz4ADAB/Gqj85l5BDgCMDEx0fA1JEnaJGOZ+SJAZr4YEXsaDYqIq4H/DFwLfHK11TO/hFzdwlL1wq6LU7MV5harRAS7B/r4hTcMM1buZ8jWRUkCmivQJoGrlx3vA16XnCLiXcCngF/NzLnWhCdJ0sZFxN8A4w1OfarZ18jM54FfjIg3AN+KiG9m5lSDcX4Jucy5+cULuy6eOjNPNZO+7tqui2PlEnuGSvT1uCApSSs1U6A9DhyMiGuAfwHuAH5z+YCIeCvw58CtmTnd8iglSdqAzHzXauciYioi9tZXz/YCF81fmflCRDwJ/Gvgmy0OdcvLTF4+O8/UbG2lbLZSa10cKtVaF8eHa62LEbYuStLFrFmgZeZiRNwNPEptF6sHMvPJiLgXOJaZR4HPAYPAN+p/eJ/LzNvbGLckSZfqKPAh4L76f7+9ckBE7ANOZeb5iLgSeCfwv2xqlAW2sFRl+nRtx8WTp3/eujgy0MeNVw0zVi4x2O8tVyVpPZr6q5mZjwCPrHju08ser/oNpSRJBXUf8HBE3AU8B7wPICImgN/OzI8AbwH+OCKS2jXZ/3Nm/j+dCrgIzs4tcmK2wtRMhZfOzpMXWhdLjA+XGB3st3VRki6BX2tJkralzDwF3NLg+WPAR+qP/y/gFzc5tELJTE6dnb+wwcfpyiIA5VIv144OMlbuZ5eti5LUMhZokiTpNeYXq0yfrhVk07NzzC9V6Ypg92AfB3bXtsIfsHVRktrCv66SJIkzc7VdF6dmK5yqty7293QxNlxivFxidKif3m5bFyWp3SzQJEnahqrV17Yunpmrty7u6OXgnkHGyiWu3Nlr66IkbTILNEmStom5xSWm69vgT5+eY6Heujgy2M+b6teT7ezznwaS1En+FZYk6TJ2urLA1GyFEzNzvHzu1dbFbvYO/3zXxR5bFyWpMCzQJEm6jFSryUtn55ierd2f7Ox8rXVxeEcv140NMjZU4gpbFyWpsCzQJEna4l5tXTwxW+FkvXWxu6vWunht/XqyHX3dnQ5TktQECzRJkrag2coCUzMVTsxWePnsPACl3m7ecMUOxsslRgb7bF2UpC3IAk2SpC2gWk1eOjPHVH2l7Fy9dfGKnX1cP15mrNzP8A5bFyVpq7NAkySpoCoLr21dXKzWWhdHB/tr15OVS5R6bV2UpMuJBZokSQXz3KlzPHvqLK+cq7Uu7ujtZt+VOxgfLjEy2E93l6tkknS5skCTJKlgzi3U2hevHy8zXi4xvLO3wxFJkjaLBZokSQVz/XiZ68c7HYUkqROa2t4pIm6NiKcj4nhE3NPgfH9EPFQ//6OIONDqQCVJkiTpcrdmgRYR3cD9wHuAG4A7I+KGFcPuAl7JzGuBzwN/1OpAJUmSJOly18wK2k3A8cx8JjPngQeBQyvGHAK+Un/8TeCWcJ9fSZIkSVqXZq5Buwp4ftnxJPD21cZk5mJEzAC7gZeWD4qIw8Dh+uFcRPx0I0FvYyOsmFM1xXlbP+ds/ZyzjXlzpwNohyeeeOKliPjnS3wZP1ONOS+rc24ac14ac15W14q5eeNG/8dmCrRGK2G5gTFk5hHgCEBEHMvMiSbeX3XO2cY4b+vnnK2fc7YxEXGs0zG0Q2aOXupr+JlqzHlZnXPTmPPSmPOyuk7PTTMtjpPA1cuO9wEvrDYmInqAYeDlVgQoSZIkSdtFMwXa48DBiLgmIvqAO4CjK8YcBT5Uf/xe4LuZ+boVNEmSJEnS6tZscaxfU3Y38CjQDTyQmU9GxL3Ascw8CnwJ+FpEHKe2cnZHE+995BLi3q6cs41x3tbPOVs/52xjnLfVOTeNOS+rc24ac14ac15W19G5CRe6JEmSJKkYmrpRtSRJkiSp/SzQJEmSJKkg2l6gRcStEfF0RByPiHsanO+PiIfq538UEQfaHVPRNTFnH4+IpyLiJxHxnYjY8H0WLhdrzdmyce+NiIwIt5WluXmLiPfXP29PRsRfbnaMRdPE7+f+iPheRPy4/jt6WyfiLJKIeCAiple792XU/Gl9Tn8SEW/b7Bg7yTzZmLlwdea8xsxpjZm3Git0bsrMtv1Q21TkH4E3AX3A3wM3rBjzb4Av1B/fATzUzpiK/tPknP0asLP++GPO2dpzVh83BPwAeAyY6HTcnf5p8rN2EPgxcGX9eE+n494Cc3YE+Fj98Q3As52Ou9M/wK8AbwN+usr524C/pnZPzXcAP+p0zJs4N+bJjc/LtsyF5rxL+sxsu5xm3rro3BQ2N7V7Be0m4HhmPpOZ88CDwKEVYw4BX6k//iZwS0Q0uvH1drHmnGXm9zLzXP3wMWr3ptvOmvmcAfwh8FmgspnBFVgz8/ZR4P7MfAUgM6c3OcaiaWbOEijXHw/z+vtGbjuZ+QMufm/MQ8BXs+Yx4IqI2Ls50XWcebIxc+HqzHmNmdMaM2+tosi5qd0F2lXA88uOJ+vPNRyTmYvADLC7zXEVWTNzttxd1Kr77WzNOYuItwJXZ+ZfbWZgBdfMZ+064LqI+LuIeCwibt206IqpmTn7DPCBiJgEHgF+Z3NC29LW+3fvcmKebMxcuDpzXmPmtMbMWxvXsdy05n3QLlGjb/hW7uvfzJjtpOn5iIgPABPAr7Y1ouK76JxFRBfweeDDmxXQFtHMZ62HWkvIzdS+nf7biLgxM3/W5tiKqpk5uxP4cmb+cUT8MrV7RN6YmdX2h7dlbec8YJ5szFy4OnNeY+a0xsxbG9exv73tXkGbBK5edryP1y+bXhgTET3UllYvttx4uWtmzoiIdwGfAm7PzLlNiq2o1pqzIeBG4PsR8Sy1PuKj2+Wi6Yto9vfz25m5kJn/BDxNLbltV83M2V3AwwCZ+UOgBIxsSnRbV1N/9y5T5snGzIWrM+c1Zk5rzLy1cR3LTe0u0B4HDkbENRHRR+3i5qMrxhwFPlR//F7gu1m/Mm+bWnPO6q0Lf04tIW2H/um1XHTOMnMmM0cy80BmHqB2rcLtmXmsM+EWRjO/n9+idiE+ETFCrT3kmU2NsliambPngFsAIuIt1BLdyU2Ncus5CnywvmPWO4CZzHyx00FtEvNkY+bC1ZnzGjOnNWbe2riO5aa2tjhm5mJE3A08Sm0XmQcy88mIuBc4lplHgS9RW0o9Tu0bwTvaGVPRNTlnnwMGgW/UrxN/LjNv71jQHdbknGmFJuftUeDXI+IpYAn4ZGae6lzUndXknH0C+IuI+F1qrRAf3gb/mL6oiPg6tZaikfo1Dn8A9AJk5heoXfNwG3AcOAf8Vmci3XzmycbMhasz5zVmTmvMvLW6Iuem2AbzL0mSJElbQttvVC1JkiRJao4FmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFcSaBVpEPBAR0xHx01XOR0T8aUQcj4ifRMTbWh+mJEmtZ46TJBVNMytoXwZuvcj59wAH6z+HgT+79LAkSdoUX8YcJ0kqkDULtMz8AfDyRYYcAr6aNY8BV0TE3lYFKElSu5jjJElF09OC17gKeH7Z8WT9uRdXDoyIw9S+gWRgYOCXrr/++ha8vSRpq3riiSdeyszRTsdxEeY4SdK6XUp+a0WBFg2ey0YDM/MIcARgYmIijx071oK3lyRtVRHxz52OYQ3mOEnSul1KfmvFLo6TwNXLjvcBL7TgdSVJ6jRznCTIrOGrAAAOkUlEQVRpU7WiQDsKfLC+09U7gJnMfF3rhyRJW5A5TpK0qdZscYyIrwM3AyMRMQn8AdALkJlfAB4BbgOOA+eA32pXsJIktZI5TpJUNGsWaJl55xrnE/i3LYtIkqRNYo6TJBVNK1ocJUmSJEktYIEmSZIkSQVhgSZJkiRJBWGBJkmSJEkFYYEmSZIkSQVhgSZJkiRJBWGBJkmSJEkFYYEmSZIkSQVhgSZJkiRJBWGBJkmSJEkFYYEmSZIkSQVhgSZJkiRJBWGBJkmSJEkFYYEmSZIkSQXRVIEWEbdGxNMRcTwi7mlwfn9EfC8ifhwRP4mI21ofqiRJrWV+kyQVzZoFWkR0A/cD7wFuAO6MiBtWDPt94OHMfCtwB/AfWx2oJEmtZH6TJBVRMytoNwHHM/OZzJwHHgQOrRiTQLn+eBh4oXUhSpLUFuY3SVLhNFOgXQU8v+x4sv7ccp8BPhARk8AjwO80eqGIOBwRxyLi2MmTJzcQriRJLdOy/AbmOElSazRToEWD53LF8Z3AlzNzH3Ab8LWIeN1rZ+aRzJzIzInR0dH1RytJUuu0LL+BOU6S1BrNFGiTwNXLjvfx+haPu4CHATLzh0AJGGlFgJIktYn5TZJUOM0UaI8DByPimojoo3aR9NEVY54DbgGIiLdQS2D2d0iSisz8JkkqnDULtMxcBO4GHgX+gdpuVk9GxL0RcXt92CeAj0bE3wNfBz6cmSvbRCRJKgzzmySpiHqaGZSZj1C7OHr5c59e9vgp4J2tDU2SpPYyv0mSiqapG1VLkiRJktrPAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCqKpAi0ibo2IpyPieETcs8qY90fEUxHxZET8ZWvDlCSp9cxvkqSi6VlrQER0A/cD7wYmgccj4mhmPrVszEHgPwDvzMxXImJPuwKWJKkVzG+SpCJqZgXtJuB4Zj6TmfPAg8ChFWM+Ctyfma8AZOZ0a8OUJKnlzG+SpMJppkC7Cnh+2fFk/bnlrgOui4i/i4jHIuLWRi8UEYcj4lhEHDt58uTGIpYkqTValt/AHCdJao1mCrRo8FyuOO4BDgI3A3cCX4yIK173P2UeycyJzJwYHR1db6ySJLVSy/IbmOMkSa3RTIE2CVy97Hgf8EKDMd/OzIXM/CfgaWoJTZKkojK/SZIKp5kC7XHgYERcExF9wB3A0RVjvgX8GkBEjFBrCXmmlYFKktRi5jdJUuGsWaBl5iJwN/Ao8A/Aw5n5ZETcGxG314c9CpyKiKeA7wGfzMxT7QpakqRLZX6TJBVRZK5st98cExMTeezYsY68tySpGCLiicyc6HQcrWaOk6Tt7VLyW1M3qpYkSZIktZ8FmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVRFMFWkTcGhFPR8TxiLjnIuPeGxEZEROtC1GSpPYwv0mSimbNAi0iuoH7gfcANwB3RsQNDcYNAf8O+FGrg5QkqdXMb5KkImpmBe0m4HhmPpOZ88CDwKEG4/4Q+CxQaWF8kiS1i/lNklQ4zRRoVwHPLzuerD93QUS8Fbg6M/+qhbFJktRO5jdJUuE0U6BFg+fywsmILuDzwCfWfKGIwxFxLCKOnTx5svkoJUlqvZblt/p4c5wk6ZI1U6BNAlcvO94HvLDseAi4Efh+RDwLvAM42uhC6sw8kpkTmTkxOjq68aglSbp0LctvYI6TJLVGMwXa48DBiLgmIvqAO4Cjr57MzJnMHMnMA5l5AHgMuD0zj7UlYkmSWsP8JkkqnDULtMxcBO4GHgX+AXg4M5+MiHsj4vZ2ByhJUjuY3yRJRdTTzKDMfAR4ZMVzn15l7M2XHpYkSe1nfpMkFU1TN6qWJEmSJLWfBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFYQFmiRJkiQVhAWaJEmSJBWEBZokSZIkFURTBVpE3BoRT0fE8Yi4p8H5j0fEUxHxk4j4TkS8sfWhSpLUWuY3SVLRrFmgRUQ3cD/wHuAG4M6IuGHFsB8DE5n5i8A3gc+2OlBJklrJ/CZJKqJmVtBuAo5n5jOZOQ88CBxaPiAzv5eZ5+qHjwH7WhumJEktZ36TJBVOMwXaVcDzy44n68+t5i7gry8lKEmSNoH5TZJUOD1NjIkGz2XDgREfACaAX13l/GHgMMD+/fubDFGSpLZoWX6rjzHHSZIuWTMraJPA1cuO9wEvrBwUEe8CPgXcnplzjV4oM49k5kRmToyOjm4kXkmSWqVl+Q3McZKk1mimQHscOBgR10REH3AHcHT5gIh4K/Dn1JLXdOvDlCSp5cxvkqTCWbNAy8xF4G7gUeAfgIcz88mIuDcibq8P+xwwCHwjIv5rRBxd5eUkSSoE85skqYiauQaNzHwEeGTFc59e9vhdLY5LkqS2M79JkoqmqRtVS5IkSZLazwJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgrCAk2SJEmSCsICTZIkSZIKwgJNkiRJkgqiqQItIm6NiKcj4nhE3NPgfH9EPFQ//6OIONDqQCVJajXzmySpaNYs0CKiG7gfeA9wA3BnRNywYthdwCuZeS3weeCPWh2oJEmtZH6TJBVRMytoNwHHM/OZzJwHHgQOrRhzCPhK/fE3gVsiIloXpiRJLWd+kyQVTjMF2lXA88uOJ+vPNRyTmYvADLC7FQFKktQm5jdJUuH0NDGm0TeFuYExRMRh4HD9cC4iftrE++vnRoCXOh3EFuS8rZ9ztn7O2ca8uYPv3bL8Bua4FvB3aP2cs41x3tbPOVu/Dee3Zgq0SeDqZcf7gBdWGTMZET3AMPDyyhfKzCPAEYCIOJaZExsJertyzjbGeVs/52z9nLONiYhjHXz7luU3MMddKuds/ZyzjXHe1s85W79LyW/NtDg+DhyMiGsiog+4Azi6YsxR4EP1x+8FvpuZDb9hlCSpIMxvkqTCWXMFLTMXI+Ju4FGgG3ggM5+MiHuBY5l5FPgS8LWIOE7tm8U72hm0JEmXyvwmSSqiZlocycxHgEdWPPfpZY8rwPvW+d5H1jleztlGOW/r55ytn3O2MR2dtzblN/DzsBHO2fo5ZxvjvK2fc7Z+G56zsFNDkiRJkoqhmWvQJEmSJEmboO0FWkTcGhFPR8TxiLinwfn+iHiofv5HEXGg3TEVXRNz9vGIeCoifhIR34mIN3YiziJZa86WjXtvRGREuBMRzc1bRLy//nl7MiL+crNjLJomfj/3R8T3IuLH9d/R2zoRZ5FExAMRMb3atvNR86f1Of1JRLxts2PcCPPbxpjj1s8ct37mt40xx61P2/JbZrbth9pF1/8IvAnoA/4euGHFmH8DfKH++A7goXbGVPSfJufs14Cd9ccfc87WnrP6uCHgB8BjwESn4+70T5OftYPAj4Er68d7Oh33FpizI8DH6o9vAJ7tdNyd/gF+BXgb8NNVzt8G/DW1e469A/hRp2Nu0WfB/LaxeTPHrXPO6uPMceuYM/PbhufNHPfa+WhLfmv3CtpNwPHMfCYz54EHgUMrxhwCvlJ//E3glohodGPQ7WLNOcvM72XmufrhY9Tu3bOdNfM5A/hD4LNAZTODK7Bm5u2jwP2Z+QpAZk5vcoxF08ycJVCuPx7m9ffV2nYy8wescu+wukPAV7PmMeCKiNi7OdFtmPltY8xx62eOWz/z28aY49apXfmt3QXaVcDzy44n6881HJOZi8AMsLvNcRVZM3O23F3UKvPtbM05i4i3Aldn5l9tZmAF18xn7Trguoj4u4h4LCJu3bToiqmZOfsM8IGImKS2O+DvbE5oW9p6/+4VgfltY8xx62eOWz/z28aY41pvQ/mtqW32L0GjbwpXbhvZzJjtpOn5iIgPABPAr7Y1ouK76JxFRBfweeDDmxXQFtHMZ62HWhvIzdS+xf7biLgxM3/W5tiKqpk5uxP4cmb+cUT8MrV7aN2YmdX2h7dlbcU8YH7bGHPc+pnj1s/8tjHmuNbbUB5o9wraJHD1suN9vH4p9MKYiOihtlx6saXCy10zc0ZEvAv4FHB7Zs5tUmxFtdacDQE3At+PiGep9QAf9SLqpn8/v52ZC5n5T8DT1BLadtXMnN0FPAyQmT8ESsDIpkS3dTX1d69gzG8bY45bP3Pc+pnfNsYc13obym/tLtAeBw5GxDUR0UftIumjK8YcBT5Uf/xe4LtZv6pum1pzzuqtDH9OLXHZM73GnGXmTGaOZOaBzDxA7ZqG2zPzWGfCLYxmfj+/Re2CfSJihFpLyDObGmWxNDNnzwG3AETEW6glr5ObGuXWcxT4YH23q3cAM5n5YqeDWoP5bWPMcetnjls/89vGmONab0P5ra0tjpm5GBF3A49S2xnmgcx8MiLuBY5l5lHgS9SWR49T+2bxjnbGVHRNztnngEHgG/XrzZ/LzNs7FnSHNTlnWqHJeXsU+PWIeApYAj6Zmac6F3VnNTlnnwD+IiJ+l1obw4e3+z/KI+Lr1NqIRurXLfwB0AuQmV+gdh3DbcBx4BzwW52JtHnmt40xx62fOW79zG8bY45bv3blt9jGcypJkiRJhdL2G1VLkiRJkppjgSZJkiRJBWGBJkmSJEkFYYEmSZIkSQVhgSZJkiRJBWGBJkmSJEkFYYEmSZIkSQVhgSZJkiRJBfH/A2myAnVfeVXBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_trace(trace);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains: 100%|██████████| 2000/2000 [00:01<00:00, 1096.82draws/s]\n",
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:1044: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  output = mkl_fft.rfftn_numpy(a, s, axes)\n",
      "E0409 10:11:36.326550 4706502080 report.py:143] There were 22 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "E0409 10:11:36.327697 4706502080 report.py:143] There were 25 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "E0409 10:11:36.328660 4706502080 report.py:143] The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    }
   ],
   "source": [
    "import pymc3 as pm3\n",
    "\n",
    "with pm3.Model():\n",
    "    mu = pm3.Normal('mu', 0, 1)\n",
    "    sd = pm3.HalfNormal('sd', 1)\n",
    "    pm3.Normal('y_0', 0, 2 * sd)\n",
    "    pm3.Normal('y_1', mu, 2 * sd)\n",
    "    pm3.sample()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
